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ABSTRACT 

We use a 24 /im selected sample containing more than 8,000 sources to study the evolution of 
star- forming galaxies in the redshift range from z = to z ~ 3. We obtain photometric redshifts for 
most of the sources in our survey using a method based on empirically-built templates spanning from 
ultraviolet to mid-infrared wavelengths. The accuracy of these redshifts is better than 10% for 80% of 
the sample. The derived redshift distribution of the sources detected by our survey peaks at around 
z = 0.6 — 1.0 (the location of the peak being affected by cosmic variance), and decays monotonically 
from z ~ 1 to z ~ 3. We have fitted infrared luminosity functions in several redshift bins in the 
range < z < 3. Our results constrain the density and/or luminosity evolution of infrared-bright 
star- forming galaxies. The typical infrared luminosity (L*) decreases by an order of magnitude from 
z ~ 2 to the present. The cosmic star formation rate (SFR) density goes as (1 + z ) 4 0±0 - 2 from z — 
to z — 0.8. From z = 0.8 to z ~ 1.2, the SFR density continues rising with a smaller slope. At 
1.2 < z < 3, the cosmic SFR density remains roughly constant. The SFR density is dominated at low 
redshift (z < 0.5) by galaxies which are not very luminous in the infrared (Xtir < 10 n Lo, where 
ixiR is the total infrared luminosity, integrated from 8 to 1000 /im). The contribution from luminous 
and ultraluminous infrared galaxies (Ltir > 10 11 £©) to the total SFR density increases steadily from 
z ~ up to z ~ 2.5, forming at least half of the newly-born stars by z ~ 1.5. Ultraluminous infrared 
galaxies (Xtir > 10 12 Lq) play a rapidly increasing role for z > 1.3. 

Subject headings: galaxies: evolution — galaxies: starburst — galaxies: photometry — galaxies: 
high-redshift — infrared: galaxies 



1. INTRODUCTION 

Infrared surveys are rapidly achieving comparable sam- 
ple sizes and areal coverage as deep ultraviolet (UV) 
and optical ones, providing an important perspective 
on galaxy evolution. Ground-based measurements plus 
the Infrared Astronomical Satellite (IRAS) have demon- 
strated that star-forming galaxies are strong infrared 
sources. Although star formation in disks is also read- 
ily detected in the UV or through optical emission-lines, 
nuclear starbursts are often heavily obscured, making 
infrared measurements essential to probe them (Ken- 
nicutt 1998). Early ground-based photometry (Rieke 
& Low 1972) and IRAS (Sanders et al. 1988) also 
revealed a population of massive galaxies in the lo- 
cal Universe with extremely high rates of star forma- 
tion (SFR > 100 A4q yr" 1 ): the ultraluminous infrared 
galaxies (ULIRGs). This violent star formation is al- 
most completely undetectable in the optical and UV part 
of the spectrum due to huge attenuation by dust. To- 
gether with lower-luminosity dust-embedded starbursts, 
this type of activity accounts for up to 20% of the local 
star formation. 

The Infrared Space Observatory (ISO) showed that 
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dust-enshrouded starbursts have undergone strong evolu- 
tion from z ~ 1 to z — (Franceschini et al. 2001). From 
results in the sub-mm with SCUBA, it appears that the 
output of ULIRGs may dominate the energy density in 
the Universe at z > 2. However, the limitations in res- 
olution and sensitivity of most of the ISO surveys have 
not allowed reliable identifications of a sufficient number 
of infrared galaxies at z < 1 to estimate robust lumi- 
nosity functions. The limitations are even more severe 
for probing the ULIRG population at higher redshifts. 
At z > 1, the 8 /im polycyclical aromatic hydrocarbons 
(PAH) band is shifted out of the longest ISOCAM band 
at 15 /im, and the longer wavelength ISOPHOT bands 
have limitations in both sensitivity and angular resolu- 
tion. Thus, we have only a first-order vision of the im- 
portance of infrared-bright galaxies at < z < 1 in the 
general picture of galaxy evolution (e.g., what percentage 
of the total star formation rate density is contributed by 
ULIRGs and what for optical/UV selected star- forming 
galaxies?). We have even less information at z = 1 — 3, 
where we believe the co-moving SFR density reaches a 
maximum (Somerville et al. 2001; Lanzetta et al. 2002), 
most of the stars in galaxies were formed (Dickinson et al. 
2003; Calura & Matteucci 2003), and dynamic structures 
(bars, disks) start to have a role in galaxy evolution (Mo 
et al. 1998). 

Spitzer's band at 24 /im encompasses PAH emission to 
z > 2, making star- forming galaxies readily detectable 
(Egami et al. 2004; Le Floc'h et al. 2004). This band 
has sensitivity an order of magnitude greater than the 
ISO 15 /im band, and 16 times as many pixels. Spitzer 
also has bands at 70 and 160 /im, with similar gains over 
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ISO to those at 24 fim. Spitzer therefore provides for 
the first time the ability to survey large fields on the 
sky to adequate depth to resolve the majority of the far- 
infrared (FIR) background and to characterize the z > 1 
population of ULIRGs and starbursts. Moreover, MIPS 6 
provides a link between the population of objects being 
discovered in the sub-mm and mm, and the UV/optical 
and near-infrared (NIR) wavelengths, giving us a key new 
tool to understand galaxy evolution. 

This work is part of a series of papers where we will 
demonstrate the ability of Spitzer to unveil galaxy evo- 
lution in the < z < 3 redshift range through its IRAC 7 
and MIPS instruments. In Bell et al. (2005) and Pa- 
povich et al. (2005) we investigate the processes govern- 
ing the evolution of star-forming galaxies from z = 
to z ~ 1. In Le Floc'h et al. (2005), we study the lu- 
minosity evolution of infrared-bright sources up to z ~ 1 
using a sample of sources detected by Spitzer in the mid- 
infrared (24 /im) , and the extensive dataset in the Chan- 
dra Deep Field South. The paper uses a combination of 
spectroscopic redshifts and photometric redshifts from 
the COMB017 8 project (Wolf et al. 2004). However, 
these redshift surveys identify few galaxies at z > 1. In 
the present paper, we extend the previously mentioned 
work to 1 < z < 3. We develop a photometric tech- 
nique based on IRAC and deep optical photometry (see 
Huang et al. 2004; Wilson et al. 2004; Lc Floc'h et al. 
2004), that is able to obtain reliable redshifts for virtu- 
ally all the galaxies detected by MIPS up to z ~ 3. With 
these redshifts, we build luminosity functions in the mid- 
infrared (MIR) in several redshift bins. This will allow us 
to study the evolution of star-forming galaxies and con- 
strain the SFR history of the Universe up to z ^ 3 using 
a homogeneously selected sample and a SFR estimator 
not affected by dust attenuation. 

This paper is organized as follows. Section 2 presents 
the observations carried out with Spitzer, and the ancil- 
lary data compiled for this study. Section 3 describes 
the technique used to estimate the redshifts from broad- 
band photometry for our sample of galaxies. Further 
details about this technique are given in Appendix A. 
The main results on the photometric redshifts are pre- 
sented in Section 4.1. The luminosity function estimation 
and fitting, constraints on the cosmic star formation rate 
density, and a discussion of the contribution of galaxies 
with different SFRs and masses to the total SFR density 
of the Universe will be presented in Sections 4.2 through 
4.5. The method for estimating the luminosity functions, 
which takes into account the photometric redshift errors, 
is described in Appendix B. Finally, the conclusions are 
outlined in Section 5. Throughout this paper, we use a 
cosmology with Ho = 70 kms _1 Mpc _1 , Om = 0.3 and 
A = 0.7. All magnitudes refer to the AB system. 

2. THE SAMPLE 

2.1. Spitzer observations 

The sample used in this paper is drawn from MIPS 
24 /urn observations of the Chandra Deep Field South 
(CDFS, a = 03' l 32 m 02 s , 5 = -27°37'24", J2000) and 

6 Multiband Imaging Photometer for Spitzer 

7 Infrared Array Camera on Spitzer. 

8 Classifying Objects by Medium-Band Observations, a spec- 
trophotometric 17- filter survey 



the Hubble Deep Field North (HDFN, a = 12 h 37 m 57 s , 
S = +62°23'14", J2000). In each field, we used the scan 
map mode to observe a rectangle of 1.5° x 0.5° in each of 
the three MIPS wavelengths (24, 70, and 160 /an). The 
overlay zone covered with all three channels is 1.0° x 
0.4°. To have the most ancillary data for each 24 /xm 
selected source, we concentrated this work in a smaller 
area around the COMB017 (Wolf et al. 2004), GOODS 9 
ACS (Giavalisco et al. 2004b), and ESO Imaging Survey 
(EIS, Arnouts et al. 2002) pointings in the CDFS case, 
and around the GOODS ACS footprint in the HDFN. In 
both fields, we also obtained IRAC data. The CDFS and 
HDFN were observed in the four IRAC channels (at 3.6, 
4.5, 5.8, 8.0 /jm) covering an area of 1.0° x 0.5° in each 
field. 

The reduction of the 24 /jm images was carried out 
with the MIPS Data Analysis Tool (Gordon et al. 2005). 
The final images had an average exposure time of ~ 
1400 s. Source detection and photometry were carried 
out with several tasks (daofind, phot, and allstar) in 
the DAOPHOT package in the Imaging Reduction and 
Analysis Facility, IRAF 10 . Sources were detected in two 
passes to recover the faintest sources, many of which are 
hidden by brighter ones. Photometry was extracted for 
all the sources (from the two passes) together to obtain 
the best possible results in crowded regions. Given the 
large point spread function (PSF) of the MIPS 24 /xm 
channel (which produces very crowded images), all mea- 
surements were made by PSF fitting. For sources of no- 
ticeable extent, the measurement aperture was set ac- 
cordingly. For the rest, a circular aperture of size ~ 15" 
was utilized. We used an aperture correction based on 
the theoretical PSF of MIPS to correct to the total flux. 
The sky estimation was carried out in two steps, first re- 
moving the large-scale variation (due to Zodiacal light) 
and then measuring the background around each source. 

IRAC images were reduced with the general Spitzer 
pipeline, and then mosaicked. The average exposure time 
of these frames is approximately 500 s. Source detection 
and photometry was carried out with SEXTRACTOR 
(Bertin & Arnouts 1996), using the same procedure as 
Huang et al. (2004). Special care was taken with the 
deblcnding of sources, given the large density of objects 
and the marked features of the IRAC PSF. Since the 
third and fourth IRAC channels are less sensitive than 
the first and second ones, we tried to obtain fluxes for 
the faintest sources in the former by carrying out the 
detection in the latter and measuring photometry in all 
bands. Photometry was performed using a small circu- 
lar aperture (3" in diameter) and an aperture correction 
was applied to get the total flux (assumed to be the one 
corresponding to a circular aperture of diameter 24.4"). 
The aperture corrections were calculated from in-flight 
PSFs. For extended sources, a circular aperture large 
enough to capture the total signal was used. 

The catalogs for the Spitzer bands in both the 
CDFS and HDFN were cut for sky regions where 
the IRAC/MIPS coverage was the deepest, and other 
UV/optical/NIR data were available (see next Section). 

9 The Great Observatories Origins Deep Survey. 

10 IRAF is distributed by the National Optical Astronomy Ob- 
servatories, which is operated by the Association of Universities 
for Research in Astronomy, Inc. (AURA) under cooperative agree- 
ment with the National Science Foundation. 
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We detected 4373 sources in the CDFS (4<r above the 
sky level) in a area of 665 arcmin 2 ; 4593 sources were 
detected in the HDFN in 517 arcmin 2 (above 3cr of the 
sky level 11 ). In Papovich et al. (2004), we estimated that 
the final catalogs are 80% complete at F 2 4 — 83 [iJy. 

2.2. Ground based optical and near-infrared photometry 

The Spitzer images were complemented with the ex- 
tensive dataset available for both the CDFS and HDFN. 
For the CDFS, we used the publicly available optical im- 
ages (UU P BVRT) released by EIS (Arnouts et al. 2002), 
the optical fluxes from COMB017 (Wolf ct al. 2004), 
the RIz frames published by the Las Campanas Infrared 
Survey (Marzke et al. 1999), the HST/ACS bviz obser- 
vations carried out by GOODS (Giavalisco et al. 2004b), 
the near-infrared JK data released by the EIS Deep Pub- 
lic Survey (EIS-DPS, Vandame et al. 2001), JHK frames 
released by GOODS (Giavalisco et al. 2004b), and the I- 
band photometry and spectroscopic redshifts released by 
the VIRMOS-VLT Deep Survey (VVDS, Le Fevre et al. 
2004). We also used UV data taken with GALEX in two 
bands at 150 nm (FUV-band) and 230 nm (NUV-band). 
For the HDFN, the Spitzer data were complemented with 
publicly available ultra-deep optical and NIR data span- 
ning from the U- to the HK s -band (UBVRIzHK Sl Ca- 
pak et al. 2004). We also used the bviz images published 
by GOODS for the central region in the HDFN. For all 
these images, source detection and photometry were car- 
ried out with SEXTRACTOR (Bertin & Arnouts 1996). 
We refer the reader to Appendix A for more detailed in- 
formation about the data compilation carried out for this 
paper. 

2.3. Redshift ancillary data 

Spectroscopic redshifts for 1599 sources have been re- 
leased by Le Fevre et al. (2004) for the CDFS (VVDS 
redshifts). In addition, COMB017 observed this field 
with up to 17 medium- and broad-band filters to obtain 
high quality photometric redshifts (Wolf et al. 2004). For 
R < 24, COMBO 17 gives redshifts for approximately 
11,000 sources. A total of 425 galaxies in our sample in 
the CDFS have VVDS spectroscopic redshifts (9% of the 
sample), and 2118 (48%) have COMB017 photometric 
redshifts. In the HDFN, several spectroscopic surveys 
have been or are being carried out. Most of the spectro- 
scopic redshifts have been compiled by the Team Keck 
Treasury Redshift Survey (TKRS, Wirth et al. 2004) and 
Cowie et al. (2004). We have also used the photometric 
redshifts found in Fcrnandcz-Soto et al. (1999). Out of 
the total number of galaxies in the HDFN survey, 601 
sources (13%) have a spectroscopic redshift. 

2.4. Merged catalogs 

Merged catalogs in all the available bands were built 
by matching the coordinates of the 24 ^m sources to 
one (the deepest) reference optical band (B band in the 
CDFS and R band in the HDFN). A 2" search radius was 

11 We used a lower detection limit in the HDFN given the ultra 
deep optical and NIR data that we had in this field. In comparison 
with the CDFS dataset, limiting magnitudes for the HDFN images 
are 0.5 — 1.0 mag fainter. This allowed us to identify fainter MIPS 
sources which were flagged as non-spurious because they were also 
detected in the optical and/or NIR. 



used. Within this search radius, multiple identifications 
(i.e., multiple sources in the optical/NIR corresponding 
to the same 24 /im detection) were found for no more 
than 7% of the sources (in the deepest ground-based im- 
ages). At this low rate, we do not expect multiple iden- 
tifications to bias our results. 

To measure the photometry, we determined the ellipti- 
cal aperture (from isophote fitting in the reference band) 
corresponding to 2.5 times the Kron radius (which con- 
tains more than 95% of the total flux of the source, ac- 
cording to Kron 1980). In all cases the apertures were 
large enough to enclose the PSF profile. This aperture 
was translated to all the other optical and NIR images. 
This procedure allowed us to obtain integrated fluxes 
for each filter in matched apertures, and to estimate the 
color properly for each source. For sources not detected 
in the reference image, we used other optical/NIR im- 
ages as the reference (if possible). In the case of the 
IRAC and MIPS bands, where the PSF is larger than 
the object images, the integrated flux was assumed to be 
that obtained from PSF fitting. 

Figure 1 shows the 24 fim flux distribution of sources in 
our sample. We also depict the 80% completeness level 
(83 /xJy). For fluxes much lower than this value, Pa- 
povich et al. (2004) showed that an increasing fraction of 
the detected sources might be spurious. We tried to iden- 
tify the real 24 /xm sources by cross-correlating the MIPS 
positions with UV-to-MIR catalogs. A simulation of the 
source density in UV-to-MIR images revealed that for a 
random position on the sky of the 24 /im image, there 
is a ~20% probability of having a counterpart in one of 
the bluer bands within the search radius. However, for a 
24 fim source with identifications in at least three other 
bands, the probability of a spurious association is almost 
negligible (less than 3%). Thus, only the 24 /im sources 
detected in three more additional bands were considered 
as real. Using this criterion, from the 4373 sources se- 
lected at 24 /im in the CDFS, 4257 (97%) were flagged as 
non-spurious detections. In the HDFN, we confirmed the 
detection for 4385 sources (96% of the total 4593 detec- 
tions). 85% of the possibly spurious sources were below 
the 83 /iJy completeness threshold. 

Of the 4257 sources in the final merged CDFS cata- 
log, 96% were detected by IRAC in at least one channel. 
The 4% remaining objects were always near very bright 
sources which interfered with the deblending algorithm. 
In the HDFN, the percentage is 89% of MIPS sources 
detected in at least one IRAC channel. 

In the optical, an average of ~ 70% of the 24 /im se- 
lected sample is detected down to B — 24.7, V = 23.8, 
R = 23.7. This percentage rises to almost 90% for 
B = 25.2, V = 24.7, and R = 24.4 (always using the 
limiting magnitudes as defined in Appendix A). These 
statistics mean that 30% of our 24 /zm selected sample 
is part of a population of IR-bright sources, starburst 
or galaxies with active galactic nuclei (AGNs), that are 
missed by UV/optical deep surveys such as COMB017. 
Only with very deep imaging (R > 25) from large tele- 
scopes (8-10 meter class) can we detect these infrared- 
bright sources. In the NIR, ~ 30% of the sources are 
detected down to J = 22.5, K = 21.8 (EIS data). In 
the UV, 30% of the sample is detected with the GALEX 
NUV channel, and 14% with the FUV filter. 

As mentioned above, our sample selection requires 
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Fig. 1. — MIPS 24 fim flux distribution of the sources detected 
in the CDFS (open histogram limited by a gray line), HDFN (filled 
gray histogram), and the total survey (open histogram limited by 
a black line). The 80% completeness flux level is shown at F24 = 
83 /xJy (Papovich et al. 2004). 



at least four identifications (at 24 /mi and three more 
bands). Most MIPS sources are detected in at least one 
IRAC channel. However, to be included in our sample 
they typically need to be also detected in at least one 
optical band 12 . Therefore, extremely red optically faint 
galaxies, which are expected to reside preferentially at 
high redshift, may not be included in the sample. How- 
ever, the number of such optically faint (i? > 25) 24 /im 
sources is at most 3-4% of the total sample (and some of 
these 24 /mi sources may be spurious) , not large enough 
to change our results significantly (e.g., the redshift dis- 
tribution presented in Section 4.1). 

3. PHOTOMETRIC REDSHIFT TECHNIQUE 

The redshift is one of the most important parameters 
to understand a distant galaxy. Although it is best de- 
termined with spectroscopy, because of the required cost 
in observing time, photometric methods are increasingly 
used. Moreover, photometric redshifts are the only ap- 
proach for very faint samples of galaxies (R > 25 or 
I > 24), given the sensitivity of the currently available 
spectrographs. There are now a number of works on the 
technique of photometric redshifts (e.g., Bolzonclla et al. 
2000; Bem'tez 2000; Collister & Lahav 2004; Wolf et al. 
2004), and on the results obtained with them (see, e.g, 
Lanzetta et al. 1996; Gwyn & Hartwick 1996; Sawicki 
et al. 1997; Rowan-Robinson 2003; Brodwin et al. 2003; 
Babbedge et al. 2004). 

These works fully develop the use of optical and NIR 
data to obtain photometric redshifts. In this wavelength 
range, there are a wide variety of spectrophotometric 
models and templates (e.g., Coleman et al. 1980; Lei- 
therer et al. 1999; Bruzual & Chariot 2003), which have 
been tested by many authors and seem to describe ac- 

12 They can also be detected only in three IRAC bands, but 
IRAC channels three and four are less sensitive than the other 
two, which causes very faint sources to be often detected only in 
the two bluer IRAC bands. 



curately the emission from all types of galaxies. Optical 
methods suffer, however, from dust attenuation, which 
can have a dramatic effect on the derived photometric 
redshifts. 

Spitzer surveys have opened another window in the 
photometric redshift possibilities. First, the galaxies de- 
tected by Spitzer might be very different in their spectral 
energy distribution (SED) properties from the typical 
UV/optical and even NIR based surveys. Second, IRAC 
observations can go much deeper than ground-based NIR 
ones in a shorter time. Therefore, one can use up to four 
more NIR and MIR bands along with the optical ones. 
Although the optical bands often include important red- 
shift indicators, the additional bands increase the con- 
fidence in the redshift determinations (Connolly et al. 
1997; Rowan- Robinson 2003). The gain in confidence 
results largely from the detection of the 1.6 /im spec- 
tral bump (see John 1988; Sawicki 2002; Lc Floc'h et al. 
2004). An issue is that models have not been developed 
in the NIR and MIR as thoroughly as in the UV/optical. 
The reason is twofold: 1) NIR spectroscopy is a relatively 
new capability (compared to optical); and 2) at these 
wavelengths, although dust attenuation is almost neg- 
ligible, dust emission starts to dominate the integrated 
spectrum of galaxies (in the continuum and in the PAH 
spectral features). Consequently, one needs detailed ra- 
diative transfer calculations to obtain models including 
both the stellar and dust contributions. Although some 
progress has been made on this topic (e.g., Corradi et al. 
1996; Gordon et al. 1997; Devriendt et al. 1999; Takagi 
et al. 2003), the set of models available in the literature 
is far from complete in the description of the emission of 
all galaxies from the UV to MIR wavelengths. 

The lack of a complete and reliable set of models de- 
scribing the SED properties of galaxies from the UV to 
the far-infrared and radio wavelengths convinced us to 
use another approach to the problem. We built empiri- 
cal broad-band SEDs for galaxies of known redshift (the 
spectroscopic sample introduced in Section 2.3). The 
resulting templates were used to fit all the galaxies in 
the entire sample. If we have a large enough number 
of previously known redshifts, and if this training set of 
galaxies is representative of the entire sample, we should 
be able to obtain reliable redshifts (see Connolly et al. 
1995). This method is similar to a neural network tech- 
nique in the sense that we use the same photometric data 
(of galaxies with known redshift) to train the photomet- 
ric redshift algorithm (see Firth et al. 2003; Collister & 
Lahav 2004) . A detailed description of our photometric 
redshift technique is given in Appendix A, jointly with a 
discussion on the reliability of the redshift estimations. 

The companion paper by Le Floc'h et al. (2005) pro- 
vides a different type of test to our redshifts. It is 
based entirely on well-tested photometric and spectro- 
scopic redshifts up to z ~ 1.1. We have compared the 
results in that paper with those reported here. As will 
be pointed out at appropriate places in the following sec- 
tions, the agreement is excellent and demonstrates that 
the outliers do not bias our results. This agreement sub- 
stantially increases our confidence in the similar results 
we obtain at z > 1.1, where there are few spectroscopic 
or previous photometric redshifts for our galaxies (see 
discussion in Appendices A and B). 
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3.1. Stars 

Given the high Galactic latitude of the fields and the 
MIR selection of the sample, stars were not expected to 
be detected in our survey in large numbers. In the final 
catalogs, only 5 sources in the CDFS and 8 in the HDFN 
were clearly identified with stars (less than 0.1% of the 
sample), based on the continuously decreasing SED in 
the NIR and MIR. This result was also checked by study- 
ing the STAR.CLASS parameter from SEXTRACTOR 
for the galaxies detected in the optical/NIR. 

3.2. AGNs 

Galaxies with an AGN are expected to be bright in the 
MIR-FIR due to the emission of the hot dust surround- 
ing the central engine (heated by x-ray and UV photons 
coming from the black hole). Indeed, IR surveys with 
ISO and Spitzer are effective in detecting AGNs (Fadda 
et al. 2002; Franceschini et al. 2002; Alonso-Herrero et al. 
2004: Rigby et al. 2004; Lacy et al. 2004). Based on the 
observed SEDs, we have tried to identify the AGNs (at 
least the most extreme cases) in our sample by selecting 
sources with monotonically rising spectra from optical to 
MIR-FIR wavelengths following a power-law and lacking 
any spectral feature (as traced by distinguishable changes 
in the slope of the SED). The MIR-FIR emission of 
these sources should be dominated by an AGN (Alonso- 
Herrero et al. 2004; Rigby et al. 2004), and they probably 
contribute non-ncgligibly to the bright end of the IR lu- 
minosity function. In addition, the photometric rcdshift 
estimation for these objects is very uncertain due to the 
lack of marked spectral features. Given that in this paper 
we are mainly interested in the evolution of star-forming 
galaxies, we carried out a first correction for the pres- 
ence of AGNs in the survey by removing from the sample 
the sources with monotonically rising spectra (Alonso- 
Herrero et al. 2005). Approximately 5% of the galaxies 
in the total sample are within this group. This percent- 
age is lower than the expected fraction of AGNs in MIR- 
FIR surveys, estimated in the range 10 — 20% by several 
authors based on local samples (Rush et al. 1993), x-ray 
identifications consistent with an active galaxy (Brandt 
et al. 2001; Hornschemeier et al. 2001; Fadda et al. 2002; 
Rigby et al. 2004), and other methods (La Franca et al. 
2004). The sources we removed are likely to be the most 
extreme cases (Type 1 obscured AGNs) within the active 
galaxy population. The existence of star formation (co- 
existing with the AGN) contributing significantly to the 
total IR luminosity of these galaxies is not completely 
ruled out, but the high dust temperatures necessary to 
obtain a power-law in the IRAC bands seem to point to a 
predominant AGN. The luminosity functions and cosmic 
SFR density calculations in the next Section were ob- 
tained from the "AGN-purged" sample. Further analysis 
of the data will be necessary to estimate the importance 
of AGNs in Spitzer surveys with higher reliability. 

4. RESULTS 

4.1. Redshift distribution 

Figure 2 shows the redshift distributions for the sam- 
ples of 24 ^m selected sources in the CDFS (in red), 
HDFN (in green), and the total (black). These dis- 
tributions result from the convolution of the real red- 
shift distribution of galaxies, the errors introduced by 



the photometric rcdshift technique, and the (flux depen- 
dent) detection curve of our survey. The cosmic variance 
between the two fields is readily apparent. We find a 
rcdshift peak around z ~ 0.7 in the CDFS, also seen 
in the COMBO 17 optical survey (Wolf et al. 2004). In 
the case of the HDFN, there seems to be a density peak 
at z ~ 0.6 and another one at z ~ 0.9, which also co- 
incides with what was found by Fernandez-Soto et al. 
(1999) for the WFPC2-HDF original field, and the spec- 
troscopic results obtained by the TKRS team. The ap- 
parent widths of these features (Az ~ 0.3) support our 
estimate of the photometric redshift errors. The curve 
for the total survey shows that the bulk of the sources 
lie at 0.5 < z < 1.0, just in the redshift range that ISO 
has probed in recent years. However, the enhanced sen- 
sitivity of MIPS in comparison with ISO has allowed us 
to detect a fainter population at z < 1.4 (see discussion 
below), and a significant number of sources at 1 < z < 3. 
In fact, almost half of the sample lies at 1 < z < 3 
(43% of the sources in both fields). This is consistent 
with the results obtained by Le Floc'h et al. (2005) using 
COMB017 and VVDS redshifts for MIPS 24 ^m sources 
detected in the optical (55% of the sources in the CDFS 
with R < 24 lying at z < 1). Note also that 75% of the 
sample is located at < z < 1.4, for which Figure A13 
and the discussion of it directly confirm the reliability of 
our photometric redshifts. 

Figure 2 also shows model predictions (cyan and 
blue histograms) for the 24 /im MIPS detections above 
F 2 a = 83 (Chary et al. 2004; Lagache et al. 2004). 
These distributions are directly comparable with the gray 
curve, built from our sample with the sources above the 
same flux cut (the shaded area depicts the differences 
in the distribution given by the photometric rcdshift er- 
rors). The general shape of the model distributions is 
roughly similar to the observations for z < 1. There 
are only small differences in the position of the redshift 
peak, which is in any case strongly affected by cosmic 
variance, as the curves for the CDFS and the HDFN 
show. However, the density of sources in each redshift 
bin below z = 1 is very different from one model to the 
other, and from the models to the data presented in this 
paper. Lagache et al. (2004) predict more sources at low 
rcdshift (z < 0.4) than Chary et al. (2004), and a less 
marked peak at z ~ 0.9. Our data lies between the mod- 
els, presenting a 15 — 25% higher number density than 
the prediction of Lagache et al. (2004) for z < 1, and a 
50% lower density than what the models by Chary et al. 
(2004) show at the peak. 

The differences between the two models and the data 
are even larger at z > 1. According to Chary et al. 
(2004), very few galaxies should be detected at high red- 
shift. On the contrary, the Lagache et al. (2004) model 
is bimodal, predicting that half of the infrared bright 
sources are at z < 1 and half at z > 1. The percentage 
of sources lying at z > 1 according to our photo- z study 
(43%) seems to be consistent with this prediction, but 
most of them are at z < 1.5, whereas the model predicts 
that most are at z > 1.4. Lagache et al. (2004) predict 
a steep decrease in the number of sources detected at 
z ~ 1.2 and a broad maximum around redshift z = 1.8 
with a width of Az ~ 1.0. This maximum is caused 
by prominent PAH features (at wavelengths from 6 to 
10 /zm) entering the MIPS 24 /<m filter as wc move to 




photo 

Fig. 2. — Observed rcdshift distribution for the 24 fim selected sources in the CDFS (red) and the HDFN (green). The combined data for 
both fields and for all flux densities are plotted in black. The predictions according to the models of Chary et al. (2004) and Lagache et al. 
(2004) for sources with 24 fim flux densities larger than 83 fijy are plotted in cyan and blue, respectively. The shaded gray area shows 
the range of results from our work (including the uncertainties in the photometric redshifts) for the 83 fijy limit (thus, this distribution is 
directly comparable with the models). 



higher redshifts. In contrast, the observed redshift dis- 
tribution shows a exponential decay with very small and 
statistically irrelevant peaks at z > 1. We do not de- 
tect the predicted minimum in source counts at z ~ 1.2, 
although this feature could be washed out due to the er- 
rors inherent to the photometric redshift technique. The 
general shape of the observed distribution is closer to the 
prediction by Chary et al. (2004), although we obtain a 
substantially higher number density at z > 1.4. 

The uncertainties in photometric redshifts are not able 
to explain the difference between the Lagache et al. 
(2004) model and our results at high z. The modest 
portion of redshift outliers expected from our estimates 
cannot be responsible for the inconsistency, either. Com- 
pared with the models, we conclude there is a lower den- 
sity of sources at high redshift or that the sources that 
we are detecting at z > 1 do not present prominent PAH 
features in the 6 < A < 10 fim wavelength range (or 
both). These sources should be very luminous (see Sec- 
tion 4.2, Figure 5, and Le Floc'h et al. 2004), and PAH 
features could be absent or hidden by a bright continuum 
or by silicate absorption, as some recent luminosity de- 
pendent models predict (see, e.g., Chary & Elbaz 2001; 
Dale & Helou 2002). Further analyses of the SEDs of 
galaxies using the three MIPS wavelengths as well as the 



ISO bands will be necessary to explore this result. 

Papovich et al. (2004), Chary et al. (2004), and Mar- 
leau et al. (2004) all presented number counts at 24 fim. 
They all found that the peak in the differential num- 
ber counts was located at a fainter flux (0.2 — 0.4 mJy) 
than predicted by the models based on ISO 15 fim ob- 
servations (roughly, at 1 mJy). In these papers, it was 
argued that the difference implies the existence of a pre- 
viously undetected population of infrared-bright galaxies 
at z ~ 1 — 3. Using the photometric redshifts derived in 
this work, we can study the contribution to the num- 
ber counts of the galaxies in different redshift bins, as 
shown in Figure 3. As we saw in the previous Figure, 
the Lagache et al. (2004) models underpredict the num- 
ber of sources at z < 1, and overpredict the number of 
galaxies above z ~ 1. Our results seem to favor a sce- 
nario where there is a strong evolution of infrared-bright 
sources from z — to z ~ 1.0, and then the evolution 
decelerates, stops or even inverts (Chary et al. 2004). We 
will come back to this issue when we present IR luminos- 
ity functions in Section 4.3. 

4.2. Infrared Luminosities and Star Formation Rates 

Models based on IRAS and ISO data on nearby galax- 
ies can be used to estimate the total infrared (TIR, inte- 
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Fig. 3. — Number counts at 24 /im built for sources in several 
redshift ranges (not corrected for completeness). Symbols of dif- 
ferent colors (joined by dotted lines of the same color for clarity) 
are used for each redshift range. The continuous lines show the 
predictions from the Lagache et al. (2004) models for each redshift 
interval (in the same color as the data points). Black filled stars 
stand for the total number counts corrected for completeness and 
presented in Papovich et al. (2004, P04 in the figure). 



grated from 8 to 1000 /tin) luminosity (see, for example, 
Sanders & Mirabel 1996; Chary & Elbaz 2001; Dale & 
Helou 2002). The most frequently used monochromatic 
fluxes to estimate TIR luminosities are at 6.7, 12, 15 /jm 
(where the highest quality ISO observations were carried 
out), and also 12, 25, 60, and 100 /tin (the IRAS bands). 
Measurements at 12 /im have been shown to be a useful 
estimator of the TIR emission for the luminosity range 
that we are dealing with (10 5 < Ltir/Lq < 10 13 , see 
Spinoglio & Malkan 1989; Spinoglio et al. 1995; Chary 
& Elbaz 2001). This conclusion is illustrated in Fig- 
ure 4, where we have plotted the relationship between 
the monochromatic luminosity at 6.7, 12, and 15 /tin 
and the TIR emission according to the models of Chary 
& Elbaz (2001). This Figure shows that the 6.7 /tm-to- 
TIR and the 15 /mi-to-TIR correlations present different 
behaviors for different luminosity ranges. The 12 /mi 
data shows the smallest scatter from a linear correlation. 
Similar results are obtained with other models found in 
the literature, such as Devriendt et al. (1999) or Dale & 
Helou (2002). 

We estimated TIR luminosities on the basis of the 
12 /mi fluxes for all the sources in the survey. This ap- 
proach also allowed us to compare with a vast number 
of papers in the literature based on IRAS 12 /mi obser- 
vations of local galaxies. We estimated the rest-frame 
monochromatic fluxes at 12 /mi by comparing the ob- 
served SEDs with models of MIR-FIR emission. There 
are several sets of these models available in the litera- 
ture (e.g., Chary & Elbaz 2001; Dale & Helou 2002). All 
of them combine an IR continuum with PAH emissions. 
The prominence of these emissions depends on the TIR 



O 



° A 

o /* 



4T 



o 6.7 /Aim 
* 12 /Aim 
a 15 fim 



8 10 12 

log(L A [L ]) 

Fig. 4. — Correlation between the monochromatic luminosity at 
6.7, 12, and 15 (im and the TIR emission (from 8 to 1000 ttm), 
according to the models by Chary & Elbaz (2001). The best linear 
fit to the 12 ^tm-to-TIR emission is also plotted (Equation 1). 



luminosity of the source. The set of templates span a 
wide range of TIR luminosities (and, thus, PAH emis- 
sion properties). However, the models are only distinct 
for wavelengths redder than ~ 3 /tm in the case of Chary 
& Elbaz (2001), and 10 /tm in the case of Dale & Helou 
(2002), although the observed SEDs in real galaxies ac- 
tually present very different shapes at bluer wavelengths 
(see Figures A12 and A16). We chose the Chary & El- 
baz (2001) template set (covering a wider range) and 
then compared these models with the observed SEDs for 
rest-frame wavelengths redder than 3 /tm. The 12 /tm 
luminosity was assumed to be the one corresponding to 
the model which best fitted the observed data. For the 
sources at highest redshift, only one point (the one for 
the 24 /zm observation) was available. In this case, we 
estimated the 12 /jm luminosity by selecting the model 
best fitting the luminosity measured by the 24 /tm chan- 
nel. 

The equations used to get the TIR emission from the 
12 /tm luminosity (Chary & Elbaz 2001), and to obtain 
SFRs from the TIR luminosity (Kennicutt 1998) are: 



log(L TIK ) = log(0.89l&#) + 1.094 x log(i 



SFR= 1.71 x 10" 10 L 



TIR 



(1) 



(2) 



where all the luminosities are in solar units, and the SFRs 
in A4q yr _1 . 

The estimation of the TIR emission from the 
monochromatic 12 /mi luminosity, and the estimation of 
the 12 /tm luminosity itself, are subject to uncertainties 
due to photometric redshift errors and the dispersion of 
the different models (within the same library and from 
one library to another; see Papovich & Bell 2002; Le 
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Floc'h et al. 2005). These uncertainties are related to 
the heterogeneous dust properties observed in galaxies 
(PAH strength and dust temperature), even for sources 
with the same bolometric luminosity ( Armus et al. 2004) . 
One more caveat in the calculation of the TIR luminosity 
is the questionable universality of the relationship given 
in Equation 1, i.e., whether the models derived from data 
for local galaxies apply to sources at higher redshifts. 
In this sense, Elbaz et al. (2005) have analyzed a sam- 
ple of sources observed with both ISOCAM at 15 fim 
and Spitzer at 24 ^m and lying at a median redshift of 
z ~ 0.7, obtaining MIR colors that require the presence 
of PAHs for more than half of their sample. Prominent 
PAHs have also been detected with IRS 13 in galaxies up 
to z ~ 3 (Houck et al. 2005). Simulations of the observed 
IRAC color-color diagrams (Sajina et al. 2005) and mod- 
els of galaxy evolution (Xu et al. 2003; Chary et al. 2004; 
Lagache et al. 2004) also seem to suggest the presence of 
PAHs in galaxies at all redshifts. However, as we move to 
higher redshifts our survey is only able to detect galax- 
ies with high TIR luminosities (LIRGs and ULIRGs, see 
Figure 5). The MIR SEDs of these sources are probably 
dominated by the continuum emission (Chary & Elbaz 
2001), and consequently the PAH importance should de- 
crease. 

Based on all these arguments, we estimate that the 
12 fim and the TIR luminosities calculated with Equa- 
tion 1 are accurate within a factor of 2-3 for individual 
galaxies. In the following discussion, we will avoid con- 
clusions that would be affected by such errors in the TIR 
luminosities. In future works, it is very desirable to deter- 
mine better MIR-FIR SEDs of galaxies at different red- 
shifts to improve the monochromatic to TIR relation. In 
fact, the uncertainty in this correlation will dominate the 
errors quoted for the luminosities of individual galaxies 
in the following sections. However, most of the results to 
come (e.g., the luminosity functions and SFR densities) 
depend on luminosities averaged over many galaxies. In 
this case, the net errors should be reduced substantially. 

Figure 5 shows the TIR luminosities of sources in our 
survey as a function of redshift (black points). Gray 
points represent sources with fluxes above the 80% com- 
pleteness limit (83 fi Jy) . A sharp detection limit is seen. 
Very few galaxies are below this sharp limit, probably 
most of them being photometric redshift outliers. This 
plot also shows that we are able to detect galaxies with 
moderate star formation (starburst galaxies with SFRs 
of a few M.q yr _1 ) up to z ~ 1. Very few ULIRGs are de- 
tected in this redshift range due to a low number density 
and/or an insufficient area coverage. Above z = 1 and 
up to z ~ 2, approximately half of the sources we are de- 
tecting present SFRs typical of infrared-luminous galax- 
ies (LIRGs with a few tens A4q yr _1 ), and half of them 
are ULIRGs (SFR>100 A^gyr" 1 ). Above z ~ 2, the 
Hyper-luminous infrared galaxy population (HyLIRGs, 
Ltir > 10 13 L Q or SFR>1700 M e yr" 1 ) starts to be de- 
tected. As expected, at z ~ 3 the detection limit reaches 
only HyLIRGs. 

The inset of Figure 5 presents the contribution of 
galaxies with different TIR luminosities to the redshift 
distribution presented in Figure 2. This Figure is not 
corrected for completeness, i.e., the distributions are af- 

13 Infrared Spectrograph for Spitzer. 



fected by the selection effects of our survey. A similar 
figure accounting for completeness effects is shown in Fig- 
ure 10. The inset of Figure 5 shows that starburst galax- 
ies (defined as the ones with Ltir < 10 11 L@) dominate 
our survey for redshifts below z ~ 0.5, and arc not de- 
tected beyond z ~ 1. The distribution for LIRGs shows a 
clear evolution from z = 0, where very few are detected, 
to z ~ 0.4, where the number of detected LIRGs starts to 
rise rapidly. At z ~ 1, all the sources in the survey are 
LIRGs. The curve for ULIRGs clarifies the statement 
in the previous paragraph about the dominance in our 
survey of this kind of sources at z > 2. 

4.3. Mid-infrared Luminosity functions 

We constructed and fitted luminosity functions at 
12 fj,m to study the evolution of the total infrared out- 
put of the population of galaxies as a function of red- 
shift, and to put constraints on the evolution of parame- 
ters such as the typical TIR luminosity of galaxies (L*). 
The luminosity functions are weakly constrained at the 
highest redshifts. Toward low luminosities, the achiev- 
able sensitivity becomes an increasingly severe limita- 
tion, while our surveyed area is inadequate to include 
rare, very high-luminosity objects. To estimate the re- 
sulting systematic errors, we have used a variety of ap- 
proaches to fit luminosity functions to the data. First, 
wc: fitted a standard form of luminosity function, allow- 
ing parametric adjustment of the density normalization, 
the slope at faint luminosities, and L* to minimize a x 2 
likelihood estimator. Second, we used a variety of func- 
tional forms determined for nearby galaxies and forced 
fits to the data at different redshifts. We were careful in 
selecting the various fitting approaches to include cases 
that would provide upper and lower limits to the TIR 
luminosity density of the galaxy population, allowing us 
to test our conclusions. 

For all the luminosity function construction, we di- 
vided the sample into redshift bins selected to provide ad- 
equate numbers of galaxies to constrain the fits well. For 
< z < 1, we used five equal intervals (Az = 0.2). We 
used four additional intervals for 1 < z < 2.6 (Az = 0.4). 
No estimations were made for z > 2.6, given the small 
range of luminosities in our sample at such high redshift. 
The estimation of the luminosity function was carried out 
using a stepwise maximum-likelihood technique (SWML, 
Efstathiou et al. 1988, see also Willmer 1997), modified 
to take into account the uncertainties in the photomet- 
ric redshifts. The procedure is described in Appendix B 
(see also Chen et al. 2003). The results were also checked 
using the V/V max method (Schmidt 1968; Huchra & Sar- 
gent 1973). 

Once we had estimated the luminosity functions with 
the modified SWML method, we carried out a variety 
of fits (see Figures 6 and 7). For the parametric fitting 
(SCHLF fit from now on), we used a Schechter (1976) 
function 14 because it includes only three free parameters 
(normalization (f>* , faint- luminosity slope a, and typical 
luminosity L*), in comparison with the four parameters 
in other parametrizations (such as a double power-law). 
Each fit at each redshift interval was independent from 
the others. 

14 The functional form is: <f>{L)dL = |1 (^Ye'T^dL 
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Fig. 5. — Selection effect on TIR luminosities of our 24 (im survey for the entire sample (black stars) and sources with flux above 83 /^Jy 
(gray stars). The divisions for LIRGs (10 11 < L Tm < 10 12 Lq), ULIRGs (10 12 < L Tm < 10 13 Lq), and HyLIRGs (L TIH > 10 13 Lq) are 
marked with horizontal lines. The inset plot shows the redshift distribution of sources with flux above 83 fijy divided into three luminosity 
ranges: starbursts (L T ir < 10 11 Lq), LIRGs (10 11 < L Tm < 10 12 Lq), and ULIRGs (L T ir > 10 12 Lq). 



The Schechter parametrization is commonly used for 
UV/optical luminosity functions. In our case, it fits the 
data points well for z > 0.2 with a minimum of free 
parameters. It is likely that the true luminosity func- 
tion is more similar to a double power-law, as many au- 
thors have shown for the local infrared galaxy population 
(Lawrence et al. 1986; Saunders et al. 1990; Rush et al. 
1993; Serjeant et al. 2001; Takeuchi et al. 2003, among 
others). Our results in the local Universe also seem to 
fit better to a double power-law, once you complement 
them with the data from other works at the bright-end 
(see Figures 6 and 7 and the discussion below). How- 
ever, within the limited luminosity range of our data for 
z > 0.2 (where we can only estimate the number den- 
sity of sources for Ltir > 10 9 L©), there is little differ- 
ence between a Schechter curve and a double power-law. 
The extra parameter required in the double power-law is 
therefore not well justified for our fits. 

We also carried out a second set of fits, using two forms 
of double power-law functions : 1) we used the local 
luminosity function derived by Rush et al. (1993); this fit 
(RUSHLF fit from now on) presents a much steeper faint- 
end slope than the local SCHLF, and a less steep slope 
at high luminosities (see Figure 6 and the discussion of 
it); and 2) we used our own derivation of the local lumi- 

15 The functional form used is: <j>(L)dL = 

(i^) 1 " 1 "" + w) dL ( Lawrcnce et al. 1986 and Rush 
ct al. 1993). 



nosity function (OWNLF from now on), that presents an 
almost flat faint-end slope (very similar to the SCHLF 
one) and is practically identical to the local RUSHLF at 
high luminosities (see the discussion below). The com- 
parison of the three different fits to the data lets us test 
directly the contribution of high luminosity galaxies to 
the overall luminosity density of the population. With 
the high luminosity shape fixed at the double power-law 
fit, the results as a function of a also let us test the 
contribution of low infrared luminosity galaxies to the 
overall output of the population. The probed values of 
a, from a = —1.7 corresponding to the local RUSHLF 
fit to a ~ —1.0 corresponding to the SCHLF or to our 
OWNLF fit, are typical for local infrared-selected galax- 
ies (Rush et al. 1993; Fang et al. 1998; Xu et al. 1998; 
Pozzi et al. 2004). This range is also very similar to the 
values found at other wavelengths for different SFR esti- 
mators and at different redshifts (see Hopkins 2004 and 
references therein, and also Cohen 2002, Poli et al. 2003, 
Bouwens et al. 2004, and Gabasch et al. 2004a,b). 

Most contemporary models of galaxy evolution assume 
that the observed change of the infrared luminosity func- 
tion with redshift can be expressed through a number 
density and/or a luminosity evolution of the local lumi- 
nosity function. This means that the shape of the lumi- 
nosity function is conserved. The evolution is normally 
parametrized with a (1 + z) n law affecting the vertical 
axis (number density), and another power-law making 
the luminosity function slide along the horizontal axis 
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Fig. 6. — Local luminosity function at 12 (im estimated by sev- 
eral authors using IRAS, ISO, and Spitzer data. Black filled stars 
show the results from this work (obtained with the entire sample of 
galaxies in the CDFS and the HDFN). Error bars are the associated 
1-er uncertainties based on Poisson statistics (Schechter 1976) and 
the propagation of photometric redshift errors (see Appendix B). 
Our results are compared with those achieved with IRAS data (for 
non-active galaxies) by Rush et al. (1993, blue open stars) and 
those obtained with ISO observations by Pozzi et al. (2004, cyan 
crossed circles). The best fit of our results to a Schechter (1976) 
function is plotted in red, the fit to two power-laws (including the 
points from Rush et al. 1993 at L12 > 10 10 Lq) is plotted in black 
(in both cases, taking into account the luminosity function errors), 
and the double power-law fit given by Rush et al. (1993) is plotted 
in blue. 



(luminosity evolution). In the double power-law set of 
fits, we allowed either a pure luminosity evolution (of 
the local luminosity function) following a (1 + z) UL law 
(L case) , or a pure number density evolution following a 
(1 + z) nD law (D case), or a combination (L + D case), 
to test the dependencies of the results on the assumed 
evolutionary behavior. The results of the different fits 
are illustrated in Figures 6 through 10. The recovered 
luminosity function parameters jointly with the derived 
luminosity densities pl 12 (for each set of fits) are given 
in Tables 1 and 2. 

Figure 6 shows the local 12 fim luminosity function 
built with all the galaxies in our CDFS and HDFN 
surveys and the best fits to it. Our results are com- 
pared with those obtained by Rush et al. (1993) using 
IRAS data and those achieved by Pozzi et al. (2004) 
using ISO observations (see also Fang et al. 1998; Xu 
et al. 1998). Given that we only detect one galaxy 
with L12 > 1O 1O L0, we assumed Rush et al. points in 
this regime to obtain the OWNLF fit. For luminosities 
10 8 % L\2 ^ 10 10 i©, our results are very close to the 
IRAS estimation, which presents slightly larger number 
densities than ISO. For L12 S 10 8 L©, our results are 
consistent with a rather flat luminosity function down to 
L12 ~ 10 6 5 Lq, and even below, in contrast to the IRAS 
results. We find a = -1.23 ± 0.06 for the SCHLF fit 
{a = -1.17±0.07forthe OWNLF fit), which is very close 



to the results obtained for the local luminosity function 
of star-forming galaxies using other SFR estimators (see, 
e.g., Serjeant et al. 2002, Perez-Gonzalez et al. 2003c, 
Wyder et al. 2005, Budavari et al. 2005, and also Fig- 
ure 8). 

Figure 7 shows the 12 /im luminosity functions for all 
the galaxies in our CDFS and HDFN surveys in the nine 
redshift bins mentioned above (including the local func- 
tion, plotted with a smaller scale than in Figure 6), and 
the fits to the data. These fits illustrate one consequence 
of our choice of the Schechter function: the contribu- 
tion at high luminosities is minimized compared with the 
more plausible double power-law fits. Thus, the SCHLF 
fits most probably underestimate the overall infrared out- 
put from high luminosity galaxies, while the RUSHLF 
and OWNLF approaches seem to provide better fits in 
this regime. A more subtle effect is that the slope to- 
ward low luminosities is not well constrained, although it 
seems to remain at a rather flat value (—1.2 > a > —1.0) 
according to the SCHLF fits and in good agreement with 
the OWNLF local value. The RUSHLF case provide the 
worst fits of the three sets, given that it assumes a very 
steep value of a. However, the small values of the slope 
may result from incompleteness in the lowest luminosity 
bin fitted, plus the poor coverage toward low luminosi- 
ties in general. Thus, it is possible that we are under- 
estimating the number density of faint sources (i.e., the 
a parameter). In the following Section, the effect of the 
large uncertainties in the faint-end slope on the estima- 
tion of the total luminosity density will be investigated 
by comparing the two extreme cases: 1) a. ~ — 1 for the 
SCHLF or OWNLF fits, either of which must provide a 
lower limit to the total output of the population of galax- 
ies in each redshift bin; and 2) a ~ —1.7 in the RUSHLF 
case, which must provide an upper limit for the lumi- 
nosity density, with an important contribution from low 
luminosity sources. 

The most robust result of the SCHLF fits is the steady 
growth of L* with increasing redshift. This is shown in 
Figure 8 (red stars), jointly with the evolution of the 
other luminosity function parameters for each set of fits. 
The normalization parameter, <f>*, grows slowly (or per- 
haps it is nearly constant, see the discussion below) to 
z ~ 0.8, and then it seems to fall. This observed decrease 
of cf>* at z > 1, which seems to be also found by other 
authors using Schechter fittings to luminosity functions 
of star-forming galaxies built at other wavelengths (cf. 
gray crosses in Figure 8), may be strongly influenced by 
the poor coverage at low luminosities, so it should be 
regarded with caution. We will come back to this issue 
later. The fitting errors in the slope parameter are large, 
and even without appealing to incompleteness, values of 
a ~ —1.2 are consistent with the data (and compatible 
with other estimations of the luminosity function of star- 
forming galaxies, cf. gray crosses in Figure 8). Since the 
Schechter function is also expected to underestimate the 
number of high-luminosity infrared galaxies, the para- 
metric fits are generally consistent with the forms of the 
luminosity function found locally. 

For the RUSHLF and OWNLF fits, the growth of L* 
with redshift is again a robust result, as shown by the 
green and blue stars in the first panel of Figure 8 16 . 

16 Note that there is an offset between the values of the L* 
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Fig. 7. — Luminosity functions at 12 fj,m for galaxies in several redshift ranges for both the CDFS and HDFN (black filled stars in all 
the panels). For the sake of completeness and clarity, we reproduce the 12 fim local luminosity function given in Figure 6 using smaller 
scales. The best fit to a Schechter (1976) function is plotted in red, the fit to two power-laws given by Rush et al. (1993) is plotted in blue, 
and our own double power-law fit is plotted in green (all with dot-dashed lines). The three fits of the local luminosity function are also 
shown in all the other panels (with the same color and line style). The continuous lines in the other panels show the best fits for each of 
the approaches described in the text: independent fit to a Schechter function (SCHLF), evolution of the local luminosity function given by 
Rush et al. (1993, RUSHLF), and evolution of the local luminosity function obtained in this work (OWNLF). Red continuous lines show 
the best SCHLF fit. Blue continuous lines show the best RUSHLF fit for each redshift interval after applying a density plus luminosity 
evolution. Green continuous lines show the best OWNLF fit for each redshift interval after applying a density plus luminosity evolution. 
Note that a pure luminosity evolution would give practically the same results (i.e., we are not able to break the degeneracy between the 
pure luminosity and the density plus luminosity evolution). 



TABLE 1 

Results of the Schechter (1976) fits (SCHLF) to the 12 fj,M luminosity functions. 





ai2 




lo g(0i 2 ) 


lo g(PL 12 ) 


Redshift range 




[L@] 


[Mpc -3 Mag -1 ] 


[L Mpc- 3 ] 


0.0 < z < 0.2 


-1.23 ±0.07 


9.61 ± 0.14 


-2.31 ±0.16 


7.38 ±0.06 


0.2 < z < 0.4 


-1.17 ±0.08 


9.81 ±0.05 


-2.18 ±0.07 


7.68 ± 0.02 


0.4 < z < 0.6 


-1.05 ±0.14 


10.05 ± 0.06 


-2.09 ±0.08 


7.98 ± 0.03 


0.6 < z < 0.8 


-1.03 ±0.12 


10.17 ±0.04 


-2.10 ±0.05 


8.08 ± 0.02 


0.8 < z < 1.0 


-1.06 ±0.46 


10.52 ± 0.11 


-2.37 ±0.17 


8.17 ±0.23 


1.0 < z < 1.4 


-1.16 ±0.15 


10.73 ± 0.05 


-2.50 ±0.08 


8.28 ±0.03 


1.4 < z < 1.8 


-1.15 ±0.34 


11.04 ±0.10 


-2.85 ±0.09 


8.24 ± 0.04 


1.8 < z < 2.2 


-1.05 ±0.55 


11.54 ±0.25 


-3.38 ±0.20 


8.17 ±0.12 


2.2 < z < 2.6 


-1.05 ±0.34 


11.86 ± 0.25 


-3.58 ±0.18 


8.29 ±0.11 
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TABLE 2 

Results of the fits to an evolved local luminosity function of the 12 /jm luminosity functions. 









RUSHLF 1 






OWNLF 2 






Evol. 3 




i°sW>i 2 ) 


l°g(P£ 12 ) 




l°g(</> 12 ) 


l°g(Pii 2 ) 


Redshift range 




[£©] 


[Mpc~ 3 Mag -1 ] 


[L Q Mpc-3] 


[£©] 


[Mpc~ 3 Mag" 1 ] 


[L Q Mpc- 3 ] 


0.0 < z < 0.2 


L 


9.90 ±0.05 


-3.10 ±0.07 


7.34 ± 0.06 


9.33 ±0.09 


-2.19 ± 0.10 


7.36 ±0.06 




L + D 


9.90 ± 0.05 


-3.10 ±0.07 


7.34 ± 0.06 


9.33 ± 0.09 


-2.19 ± 0.10 


7.36 ± 0.06 


0.2 < z < 0.4 


L 


10.24 ±0.07 


-3.10 ±0.07 


7.66 ± 0.03 


9.59 ±0.10 


-2.19 ± 0.10 


7.57 ±0.02 




L + D 


10.09 ±0.06 


-2.95 ±0.09 


7.65 ± 0.02 


9.55 ±0.10 


-2.14 ±0.10 


7.62 ± 0.02 


0.4 < z < 0.6 


L 


10.54 ±0.12 


-3.10 ±0.07 


8.04 ± 0.05 


9.81 ±0.12 


-2.19 ± 0.10 


7.87 ± 0.03 




L + D 


10.25 ±0.09 


-2.82 ±0.12 


8.03 ± 0.04 


9.73 ±0.11 


-2.11 ± 0.11 


7.85 ± 0.03 


0.6 < z < 0.8 


L 


10.79 ±0.16 


-3.10 ±0.07 


8.21 ±0.06 


10.01 ±0.14 


-2.19 ± 0.10 


8.02 ± 0.03 




L + D 


10.39 ±0.12 


-2.70 ±0.15 


8.21 ±0.04 


9.90 ±0.13 


-2.07 ±0.13 


8.04 ± 0.02 


0.8 < z < 1.0 


L 


10.85 ±0.20 


-3.10 ±0.07 


8.29 ± 0.24 


10.11 ±0.17 


-2.19 ±0.10 


8.15 ±0.24 




L + D 


10.61 ±0.15 


-2.84 ±0.18 


8.31 ±0.24 


10.16 ±0.15 


-2.24 ±0.14 


8.15 ±0.23 


1.0 < z < 1.4 


L 


10.93 ±0.25 


-3.10 ±0.07 


8.36 ±0.10 


10.23 ±0.20 


-2.19 ± 0.10 


8.26 ± 0.06 




L + D 


10.74 ±0.19 


-2.86 ±0.22 


8.42 ± 0.07 


10.35 ±0.17 


-2.37 ±0.16 


8.20 ± 0.04 


1.4 < z < 1.8 


L 


11.02 ±0.30 


-3.10 ±0.07 


8.46 ±0.12 


10.38 ±0.24 


-2.19 ± 0.10 


8.41 ±0.07 




L + D 


11.28 ±0.23 


-3.47 ±0.27 


8.35 ± 0.08 


10.82 ±0.21 


-2.86 ±0.18 


8.18 ±0.04 


1.8 < z < 2.2 


L 


11.17 ±0.35 


-3.10 ±0.07 


8.60 ±0.18 


10.55 ±0.28 


-2.19 ± 0.10 


8.58 ±0.14 




L + D 


12.08 ±0.27 


-4.23 ±0.31 


8.38 ±0.15 


11.33 ±0.24 


-3.36 ± 0.20 


8.19 ±0.12 


2.2 < z < 2.6 


L 


11.42 ±0.39 


-3.10 ±0.07 


8.86 ±0.18 


10.80 ±0.31 


-2.19 ± 0.10 


8.83 ±0.13 




L + D 


12.21 ±0.30 


-4.18 ±0.35 


8.56 ±0.15 


11.63 ±0.26 


-3.46 ±0.22 


8.39 ±0.11 



Note. — 1 For the RUSHLF fits, the slopes at low and high luminosities arc fixed: a 12 = -1.70 ± 0.02 and 12 = -3.60 ± 0.09. 2 For the 
OWNLF fits, the slopes are also fixed: Q12 — —1.17 ± 0.07 and /3i2 — —2.97 ± 0.16. 3 The type of evolution can be: luminosity evolution (L), 
or a combined luminosity plus number density evolution (L + D). 
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This result is obtained for both pure luminosity evo- 
lution (open stars) and combined luminosity plus den- 
sity evolution (filled stars), although the rate of change 
varies from one type of evolution to the other (and from 
these to the SCHLF approach) . The rate of change from 
z = to z ~ 1 seems to follow a (1 + z) nL law, where 
til = 3 — 5. In the region of overlap (z < 1.1), this 
behavior agrees closely with the results of Le Floc'h ct 
al. (2005). The uniform pure density evolution seems 
to be ruled out by our data, given that it fails to de- 
scribe the luminosity function points as early asz~ 0.4. 
Indeed, the probability of exceeding by chance the x 2 
value obtained for the fit with a pure density evolution is 
7% for the RUSHLF fits (5% for OWNLF). In contrast, 
the probability of exceeding by chance the \ 2 value ob- 
tained for the fit with pure L or L + D evolution is 75% 
(80%) and 90% (93%), respectively. This means that wc 
are not able to confidently break the degeneracy between 
these two scenarios, although the fits for the combined 
L + D evolution are slightly better. That is, the luminos- 
ity function data points can be reproduced with a strong 
L evolution (til ~ 4 — 5), or with a weaker luminosity 
evolution (til — 2.6 — 3.1) combined with a relatively 
smaller density evolution [no = 0.5 — 2.0). When we 
used Schechter functions, which also allowed changes in 
the faint end slope, the combined L + D evolution is also 
slightly favored. 

Figure 8 also shows the best linear fits to the evolution 
of the luminosity function parameters at < z < 0.8. In 
the case of the Schechter fitting, the number density of 
sources, as parameterized by <fi* , evolves as (1 + z )°- 9±0 - 6 
in this redshift range. The typical infrared luminos- 
ity (L*) evolves as (1 + z ) 3 - 1±6 - 5 . In the OWNLF fits, 
the L + D evolution predicts L* oc (1 + z ) 3 0±0 - 3 and 
<p* oc (1 + 2)°- 6±0 - 2 (green dashed lines in the first two 
panels) . In comparison, when we fit the data with a pure 
L evolution in the OWNLF case, 4>* remains obviously 
constant, and L* behaves similarly to the Schechter func- 
tion fits, evolving as (1 + z ) 3 - 6±a - 3 (green dotted lines 
in the first two panels). As we mentioned before, the 
RUSHLF fits are always considerably worse than the 
OWNLF and SCHLF ones, but they predict similar evo- 
lution laws: L* oc (1 + z ) 2 -6±i.i an d 0* oc (1 + z) 2A±0M 
for the L + D case, and L* oc (1 + z ) 4 - 7±0 - 3 for the pure 
L evolution. Combining the results from the three sets 
of fits, the most probable values of the exponents of the 
evolution laws are: ni = 3.0 ± 0.3 and no = 1.0 ± 0.3. 

Above z ~ 0.8, the evolution degeneracy is more se- 
vere. The Schechter fits seem to favor a scenario where 
<p* decreases steadily as (1 + z )- 51 +°- 8 j an d L* contin- 
ues rising as (1 + z )4.8±o.8^ j e ^ a f ew g a i ax i es w ith very 
violent star formation dominate the TIR luminosity den- 
sity. However, our data could also be reproduced with a 
slower (or even null) decrease of </>* up to z ~ 3 and an 
also slower increase of L* (see open stars in the left two 
panels of Figure 8). 

4.4. Cosmic star formation rate density 

In the previous Section, we presented two ways of fit- 
ting the luminosity functions: one using a Schechter 
function to fit the data independently at each individ- 
ual redshift range, and the other assuming a constant 
shape of the luminosity function (the local shape, which 
is the best constrained) and evolving it with redshift 



in density or/ and luminosity. The former technique al- 
lows a change in the shape of the luminosity function, 
mostly in the faint end slope, while the latter proce- 
dure fixes this slope. It also assumes a less rapid fall 
at high-luminosities, which seems to be the case in the 
local Universe. In this Section, we integrate the SCHLF, 
OWNLF, and RUSHLF fits of the luminosity functions 
to get different (possibly biased) estimates of the total 

12 /im luminosity density of the Universe at < z < 3, 
the TIR luminosity density, and the SFR density. The 
TIR luminosities have been estimated using Equation 1 . 
These luminosity densities will be translated to cosmic 
SFR densities in Figure 9 using Equation 2. At the end 
of this Section, by comparing the results obtained with 
the three sets of fits, we will discuss how changes in the 
faint end and bright end slopes affect the estimations of 
the luminosity and SFR densities. 

Figure 9 shows the evolution of the SFR density of the 
Universe as a function of redshift (red, green, and blue 
stars referring to the integration of the SCHLF, OWNLF, 
and RUSHLF fits, respectively) in a Lilly-Madau dia- 
gram (Lilly et al. 1995; Madau et al. 1996). Our survey 
reproduces the rapid increase in psfr from z = to 
z ~ 1.4 observed by many previous works. Our esti- 
mations follow a (1 + 2;) 4 0±0 - 2 law up to z = 0.8 (i.e., 

13 = 4.0 ±0.2), and a lower slope {(3 - 3.4) up to z ~ 1.4. 
For z < 1, we obtain: 



log (psfr) = (-1.87 ± 0.04) + (3.98 ± 0.22) x log (1 + z) 

(3) 

This result is consistent with that of Hopkins (2004, 
see also Hogg 2002), who used all the SFR density es- 
timations plotted in Figure 9 (obtained with different 
SFR tracers), although our slope is higher (and our lo- 
cal density is slightly smaller): Hopkins (2004) gives 
(3 = 3.10 ± 0.25 applying a simple obscuration correction 
for the SFR density estimations, and (3 — 3.29 ± 0.26 
for a luminosity-dependent obscuration correction. The 
analysis of Hopkins (2004) is partly based on UV surveys, 
which tend to obtain significantly less steep values of the 
SFR density evolution slope ((3 = 2 — 2.5, sec Schimi- 
novich et al. 2005; Baldry et al. 2005). However, other 
works based on UV and optical surveys (using different 
emission-lines) find larger values, closer to our estimation 
(e.g., Lilly et al. 1996; Tresse et al. 2002; Hippelein et al. 

2003) . Moreover, models based on IR and sub-millimeter 
models also predict an evolution with an exponent close 
to (3 ~ 4 (Blain ct al. 1999; Xu ct al. 2003; Lagache et al. 

2004) . 

The scatter in the results for the slope of the evolution 
of Psfr suggests that the extinction properties of the 
galaxies dominating the total SFR density are evolving 
with redshift. For example, our estimation of the SFR 
density at z ~ 0.1 (which is consistent with other esti- 
mations based on radio observations, e.g., Condon 1989; 
Sadler et al. 2002; Serjeant et al. 2002) is ~40% lower 
than the most recent results achieved by Ha (Perez- 
Gonzalez et al. 2003c) or UV/optical surveys (Glaze- 
brook et al. 2003; Wyder et al. 2005; Schiminovich et al. 
2005; Martin et al. 2005) in the local Universe. This may 
indicate that the star formation in this redshift regime 
is dominated by galaxies with not very extincted bursts, 
where the dust emission only traces a small part of the 
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Fig. 8. — Evolution of the 12 /im luminosity function parameters from z = to z ~ 3. For all panels, the results from the SCHLF, 
RUSHLF, and OWNLF fits are plotted in red, blue, and green, respectively. Open stars refer to the results for a pure luminosity (L) 
evolution, while filled stars refer to a combined luminosity plus density evolution (L + D). The best fits for the evolution of L* and 0* at 
z < 0.8 are plotted with dotted and dashed lines for L and L + D evolutions, respectively, with the same color code as the points for the 
three sets of fits. Gray symbols show some comparison values for the three parameters (for a Schechter parametrization) extracted from 
the literature and based on fits of luminosity functions built with samples of star-forming galaxies selected with different SFR estimators 
(Gallego et al. 1995; Connolly et al. 1997; Tresse & Maddox 1998; Cowie et al. 1999; Steidel et al. 1999; Yan et al. 1999; Machalski & 
Godlowski 2000; Moorwood et al. 2000; Hopkins et al. 2000; Sullivan et al. 2000; Wilson et al. 2002; Gallego et al. 2002; Sadler et al. 2002; 
Serjeant et al. 2002; Tresse et al. 2002; Teplitz et al. 2003; Perez-Gonzalez et al. 2003c; all these estimations were compiled by Hopkins 2004, 
from which we extracted the SFR calibrations to convert all the L* values based on different estimators to L^ 2 )- The evolution laws plotted 
in the Figure for a L + D scenario are: L* evolves as (1 + z)3.1±0.5 and ^* evo i ves as (i + ^,)0.9±0.6 for the SCHLF fitting (red dashed 
lines); L* tx (1 + z)2-6±l.l an d 0* oc (1 + z ) 2 - 1±0 - 6 for RUSHLF (blue dashed lines); and L* oc (1 + z) 3 - 0±0 - 3 and 0* oc (l + z) a6±0 - 2 for the 
OWNLF case (green dashed lines). For a pure L evolution: L* oc (1 + 2:) 4 7±0 - 3 for RUSHLF (blue dotted line); and L* oc (1 + z )3 6±0.3 
for OWNLF (green dotted line). 



total SFR of each galaxy (i.e., many photons from the 
newly-born stars do not interact with the dust, but they 
escape through the UV or emission-lines). This effect 
would be reasonable in galaxies with not very intense 
star formation (SFR < 5.Moyr _1 ), given that the most 
violent star-forming galaxies show the highest dust atten- 
uations (Hopkins et al. 2001; Sullivan et al. 2001; Perez- 
Gonzalez et al. 2003c). Those low intensity star-forming 
galaxies contribute importantly to the total SFR density 
in the local Universe (see Figure 16 in Perez-Gonzalez 
et al. 2003b). As we move to higher redshifts, the star 
formation in galaxies starts to be dominated by intense 
dust-enshrouded bursts, whose SFR is better traced by 
the TIR emission (Cardiel et al. 2003). In this scenario, 
the evolution of the extinction properties also seems to be 
consistent with the increasing contribution of IR-bright 
galaxies with dust enshrouded bursts to the cosmic SFR 
density (Chary & Elbaz 2001; see also Figure 10). 

After the increase from z = to z ~ 1.4, we find a 
roughly flat behavior of Psfr up to z ~ 3, very similar 
to what some models predict (e.g., Lagache et al. 2004, 
shown in the figure, or Chary & Elbaz 2001), and consis- 
tent with the results from most UV/optical surveys (see 
Hopkins 2004, and references therein). 

There are three issues affecting our estimations in the 
Lilly-Madau diagram: how the fitting procedures affect 
the results, the translation of 12 /im luminosity to bolo- 
metric infrared luminosity, and of bolometric infrared lu- 
minosity to true bolometric luminosity and star forma- 
tion rates. 

For the first of the issues, by integrating under the lu- 
minosity functions of a given shape, we found that the 
TIR output of the galaxy population was virtually inde- 



pendent of whether the evolution was in luminosity, den- 
sity, or a mixture (within the uncertainties; see Table 2). 
By varying the shape of the fitting curve, we also found 
that the total output of the galaxy population was not 
strongly dependent on the shape of the luminosity func- 
tion above L* . This is shown with green stars in Figure 9. 
These points were calculated using the OWNLF fits, i.e, 
an almost flat behavior at the faint end and a less steep 
behavior than the Schechter function at the bright end. 
It appears that the luminosity function is too steep in the 
high-luminosity region for plausible variations to change 
the integral of the luminosity function significantly (aver- 
age change smaller than 20%). Given that AGNs should 
predominantly populate the bright end of the luminosity 
function, it seems probable that they do not affect our 
results by a large factor, either 17 . 

For a fixed behavior toward high luminosity, we found a 
change by a factor of two as a was changed from —1.7 (as 
in the RUSHLF case and some UV surveys) to approx- 
imately -1.0 (as the SCHLF and OWNLF estimations 
predict). The range of estimates is shown in Figure 9 
with the shaded area. Therefore, the uncertainties in 
the TIR luminosity density include a significant contri- 
bution from the lack of knowledge of the low luminosity 
galaxy population. At low redshift, a seems to be close 
to the fiat value: -L3<a<-1.0 (Serjeant et al. 2002, 
Perez-Gonzalez et al. 2003c, Wyder et al. 2005, Budavari 
et al. 2005). The most recent estimations of a based on 
the deepest observations of high-redshift galaxies emit- 
ting strongly in the UV (Gabasch et al. 2004a) suggest 

17 For example, the AGNs removed from our sample as discussed 
in Section 3.2. 
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Z 

Fig. 9. — The Lilly-Madau diagram (Lilly et al. 1995; Madau et al. 1996): evolution of the SFR density of the Universe with redshift. The 
estimations based on the SCHLF, RUSHLF, and OWNLF fits are plotted with red, blue, and green stars, respectively (see text for details). 
The shaded area delimits the zone between the two extreme SFR density estimations for each redshift. The heavy error bar at z ~ 0.1 
shows the uncertainty in the transformation from the monochromatic 12 (im luminosity to the TIR emission, as shown in Equation 1 (sec 
the text for a discussion of this error). This error is common for all our points, but it is only given in the first one for clarity. Vertical 
segments for each point show the uncertainty related to the integration of the SCHLF luminosity function (comparable to the OWNLF 
and RUSHLF cases). The horizontal lines show the range of rcdshifts used in each bin. The curves show two typical models: one with 
a decay from z ~ 1 (Xu et al. 2003, SI model), and another with a constant SFR density at high redshift (Lagache et al. 2004). The 
colored points (shown with error bars) are extracted from different sources in the literature, normalized to the same cosmology by Hopkins 
(2004). Red symbols are estimations based on Ha or H/3 measurements (Gallego et al. 1995; Pettini et al. 1998; Tresse & Maddox 1998; 
Glazebrook et al. 1999; Yan et al. 1999; Moorwood et al. 2000; Hopkins et al. 2000; Sullivan et al. 2000; Tresse et al. 2002; Perez-Gonzalez 
et al. 2003c). Green symbols stands for [0//]A3737 estimations (Hammer et al. 1997; Hogg et al. 1998; Gallego et al. 2002; Teplitz et al. 
2003). UV-based data points are plotted in blue (Lilly et al. 1996; Connolly et al. 1997; Treyer et al. 1998; Cowie et al. 1999; Steidel et al. 
1999; Sullivan et al. 2000; Massarotti et al. 2001; Wilson et al. 2002; Giavalisco et al. 2004a). Cyan estimations are based on mid-infrared 
data (Flores et al. 1999). Magenta points are based on sub-mm and radio observations (Condon 1989; Hughes et al. 1998; Barger et al. 
2000; Machalski & Godlowski 2000; Haarsma et al. 2000; Sadler et al. 2002; Serjeant et al. 2002; Condon et al. 2002). The yellow point is 
based on X-ray data (Georgakakis et al. 2003). 



also an almost flat value with a marginal indication of 
evolution of the slope with redshift (to shallower values). 
This means that the SFR densities should be closer to 
the SCHLF or OWNLF values than to the upper lim- 
its of the RUSHLF case (quoted with the blue stars in 
Figure 9). 

The very weak dependence of the integral of the lu- 
minosity function on either very high or very low (for 
a ~ — 1) luminosity galaxies also means that photomet- 
ric redshift outliers have little effect on this integral. In 
addition, the integral of the luminosity function is always 
a more robust calculation than the individual parameters 
of the fitting function. The Monte Carlo simulation de- 
scribed in Appendix B confirms these statements. 

The second of the issues affecting our SFR den- 
sity estimates is related to the uncertainties in the 
monochromatic-to-TIR emission relationship. This error 
is roughly a factor of two for individual galaxies (based 
on Equation 1, extracted from Chary & Elbaz 2001). We 



show this range on the first point (at z ~ 0.1) in Figure 9. 
Fortunately, these uncertainties can be reduced by future 
work including the longer wavelength MIPS bands and 
ground-based data in the sub-millimeter. Moreover, for 
calculating the TIR luminosity density, we average the 
luminosities of galaxies in luminosity bins, and integrate 
them to all the possible values. In this averaging and 
integration procedure, it is probable that the final un- 
certainties diminish. 

The third issue in the SFR density estimation appears 
to be less of a problem as we move to higher redshifts. 
The unaccounted contribution of ultraviolet luminosity 
is probably lower than 40% of the infrared contribution 
for all redshifts. In this regard, Bell et al. (2005) quote 
that the UV luminosity density is ~ 7 times smaller than 
the IR density at z ~ 0.7. In Burgarella et al. (2005), an 
average UV attenuation of Afuv = 2.7 for a sample of 
ULIRGs at z < 1.6 is found. Attenuations of the same 
order (factors of 5 — 10) have also been found by other au- 
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Fig. 10. — Relative contribution of starbursts (Ltir, < 10 11 Lq), 
LIRGs (10 11 < Ltir < 10 12 Lq), and ULIRGs (Ltir > 10 12 L ) 
to the total SFR density of the Universe as a function of rcdshift. 
The two extreme cases of the faint-end slope value are separated: 

a ~ -1.7 (RUSHLF fits, upper panel) and a 1.0 (SCHLF 

case, lower panel). The error bars show the uncertainties on the 
integration in each luminosity range, considering the errors of the 
individual luminosity function parameters. 



thors (Madau et al. 1998; Steidel et al. 1999; Adelberger 
& Steidel 2000; Massarotti et al. 2001; Schiminovich et al. 
2005, among others). This means that the star formation 
traced by the UV alone (without extinction correction, 
i.e., the star formation which the IR cannot trace because 
it did not heat the dust) is 5 — 10 times smaller than the 
star formation traced by the IR. If there is an evolution 
in the extinction properties of galaxies, as we previously 
discussed, these uncertainties will yield an offset of the 
total SFR density which depends on the rcdshift. If the 
evolution is not present, there should be a systematic 
offset of the estimates (i.e., all the points in the Lilly- 
Madau plot would move up approximately by the same 
quantity). A comparison of the extinction properties of 
the galaxies selected in UV/optical and IR surveys will 
be necessary to address this issue. 

4.5. Contribution of galaxies with different TIR 
luminosities and masses to the total SFR density of 
the Universe 

Figure 10 shows the contribution of galaxies of different 
TIR luminosities to the integrated SFR (or TIR luminos- 
ity) density of the Universe as a function of redshift. We 
have produced two plots to account for the two extreme 
cases of the faint-end slope value: a ~ —1.0 (SCHLF 

case) and a 1.70 (RUSHLF fits). In the almost flat 

luminosity function scenario, there is a dominant but de- 
creasing contribution of normal and starburst galaxies 
with faint infrared luminosities (Ltir < 10 11 L@) to the 



total SFR density up to at least z ~ 0.8. At this red- 
shift, LIRGs already form approximately half of the total 
amount of newly-born stars. These results are consistent 
with the ones achieved by Le Floc'h et al. (2005) using 
CDFS data and photometric redshifts from COMBO 17. 
The evolution of LIRGs was expected by ISO-based mod- 
els, such as Chary & Elbaz (2001); Chary et al. (2004) 
and Xu et al. (2003), which predicted a ~ 70% contribu- 
tion to the total luminosity density for z > 0.5. Our es- 
timation is somewhat below this value. The evolution of 
LIRGs decelerates at z ~ 0.9, remaining approximately 
at the 50 — 60% level up to z ~ 1.5, while starbursts 
continue their decline, and ULIRGs start to contribute 
significantly to the total luminosity density. By z ~ 2, 
ULIRGs already form more than 70% of the newly-born 
stars in the Universe, and they completely dominate the 
luminosity density at z ~ 3. 

Figure 10 also shows the contributions to the total SFR 
density obtained with the RUSHLF case, i.e., consider- 
ing a rather steep luminosity function for all redshifts 
(upper panel). As shown in Figure 9, this gives an up- 
per value for the SFR density of the Universe, with an 
important contribution from galaxies with modest star 
formation. Indeed, the upper panel of Figure 10 shows 
that starbursts still show a dominant but decreasing con- 
tribution to the total SFR density at z < 0.6, but their 
contribution is not negligible at z ~ 2.0 (they still form 
roughly 30% of the total amount of stars at that redshift). 
The contribution of LIRGs rises slowly up to ~ 30% at 
z ~ 1.5, and then stays approximately constant up to 
z ~ 3.0. Starting at z ~ 1, ULIRGs start to contribute 
non-negligibly to the total SFR density, and reach ~30% 
of the total SFR density at z ~ 2. Note that the contri- 
bution from starbursts to the total SFR density is larger 
in the flat slope case than in the steep case for z < 0.3 
(and, consequently, the contribution of LIRGs to the to- 
tal density is smaller in the a ~ —1.0 case). This effect 
is directly related to the uncertainties in the fits (in L* 
and a) at low redshift. 

Figure 10 demonstrates the shift of the star forma- 
tion density to LIRGs and ULIRGs (i.e., to IR-bright 
galaxies with very violent dust enshrouded bursts of star 
formation) as we move from z = to z > 1. LIRGs 
and ULIRGs tend to be the most massive star-forming 
galaxies at z < 1, with masses M* ~ 10 11 Mq (see, e.g., 
Rigopoulou et al. 2002; Tacconi et al. 2002; Franceschini 
et al. 2003). It is interesting to study the connection 
between dust enshrouded star formation and the stellar 
mass of each galaxy. We analyze the relationship between 
these two parameters in Figure 11. 

We estimated stellar masses for all our galaxies using 
K-b&nd luminosities calculated by interpolation among 
the templates utilized to get the photometric redshifts. 
By using the IRAC photometry, we could probe the 
rest- frame if -band up to z ~ 3. For local galaxies, a 
number of authors have demonstrated how accurate stel- 
lar masses can be determined from if-band luminosities 
combined with optical colors to constrain the star for- 
mation history (Bell & de Jong 2001; Bell et al. 2003; 
Kauffmann et al. 2003). They show that the mass- to- 
light ratios for if-band luminosities should not change 
more than a factor of 2 to 3 across a wide range of star 
formation histories, in comparison with a factor of >10 
for optical mass-to-light ratios. The effects of extinction 



17 



are also negligible in the near infrared. 

At higher redshifts, due to galaxy evolution, the mass- 
to-light ratio should decrease (Drory et al. 2001, 2004; 
Fontana et al. 2004). We therefore based our mass es- 
timates on the redshift-dependent relationships for late- 
type galaxies found in Fontana et al. (2004). At low red- 
shift, this procedure agrees well (with less than a 10% 
scatter) with, for example, that of Bell et al. (2003). 
However, for 0.2 < z < 3, the values are down to 2 
times lower than the local relationship found by Bell et al. 
(2003). A lower limit to the masses can also be deter- 
mined by starburst modeling, based on the galaxy lumi- 
nosities. The models of M82 (Rieke et al. 1993; McLeod 
et al. 1993; Forster Schreiber et al. 2003) suggest that the 
star formation in the infrared-luminous galaxies would 
produce a stellar population only slightly (no more than 
a factor of two) lower in mass than the value determined 
by the method of Fontana et al. (2004). We conclude 
that our mass estimates are accurate to a factor of 2 to 
3, which is sufficient for the following qualitative analy- 
sis of the connection between star formation and stellar 
mass as a function of redshift. More robust estimates of 
the stellar masses (leading to a more detailed study of the 
the SFR-mass connection) should rely on the modeling 
of the stellar populations in each galaxy (see, e.g., Pa- 
povich et al. 2001; Kauffmann et al. 2003; Perez-Gonzalez 
et al. 2003a, b), a topic that will be addressed in future 
works. 

Figure 11 shows the relationship between the specific 
SFRs (SFR per stellar mass unit) and the total stellar 
masses for all the galaxies in our survey, divided into 
redshift bins. As we move to higher redshifts, there is a 
trend for the most active star-forming galaxies to be more 
massive: the cloud of points for each redshift range shifts 
to larger stellar masses and specific SFRs as we move to 
higher redshifts. The trend is in part due to selection 
effects, i.e., we are only detecting the IR brightest objects 
at the furthest distances, and those objects should also 
be bright in the optical and NIR, thus presenting large 
stellar masses. However, it is worth pointing out that 
our survey is able to detect the galaxies that dominate 
the SFR density up to z <~ 3, and those galaxies are more 
and more massive. 

This result seems to support the theory of "downsiz- 
ing" (Cowie et al. 1996), for which there is an increasing 
amount of observational evidence (Heavens et al. 2004; 
Glazebrook et al. 2004; Juneau et al. 2005; Bauer et al. 
2005). In this theory, the most massive galaxies form 
first in very violent episodes of star formation, while the 
formation of less massive systems continues as we move 
to lower redshifts. Juneau et al. (2005) argue that the 
SFR in the most massive galaxies (A4* > 10 108 M.q) was 
much larger at z ~ 2 than in the local Universe. This is 
supported by Figure 11 because most of the galaxies with 
A4* > to 10 - 6 - 10 - 8 M Q and high specific SFRs are placed 
at z > 1.5. That mass value is, in fact, a rather sharp cut- 
off for the entire sample, which should have its origin in 
the steep fall of the stellar mass function (see Drory et al. 
2004; Fontana et al. 2004; Feulner et al. 2005), but may 
also be due to the most massive systems having formed 
the bulk of their stars in an epoch before z = 1.5 — 2. 5 in 
very violent episodes of star formation, probably present- 
ing very high IR luminosities. Beyond that point, their 
specific SFRs decrease considerably, gradually disappear- 



ing from our plot. This scenario would also be supported 
by the results about the population of high-mass galaxies 
at z = 1—3, that require ULIRGs at z > 3 (see, e.g., Mc- 
Carthy et al. 2004; Labbe et al. 2005). For intermediate- 
mass galaxies (10 10 < M« < 1O 1O 6 A4 ) with high 
SFRs, the observed redshift distribution peaks in the 
range 1.0 < z < 2.5, also in agreement with Juneau et al. 
(2005). The dominant sources in our survey at redshifts 
z < 1, just when the cosmic SFR density starts its de- 
cline, present M* < 10 10 M Q and SFR/A4, < SGyr" 1 . 

5. CONCLUSIONS 

We have cross-correlated the sources detected by MIPS 
at 24 /jm in the Chandra Deep Field South and Hubble 
Deep Field North with ultraviolet, optical, near-infrared, 
and mid-infrared (IRAC) catalogs. 

Using this multiwavclcngth dataset, we have estimated 
the photometric redshifts of all the sources in our sample. 
The technique used to estimate these redshifts is based on 
empirically-built templates obtained from sources with 
know spectroscopic redshifts. The accuracy of these red- 
shifts is better than 10% for 80% of the sample. As a 
test of our conclusions based on this redshift estimation, 
we show that our results for < z < 1 closely agree with 
those in a companion paper by Le Floc'h et al. (2005). 
The derived redshift distribution of the sources detected 
by our survey (for fluxes F24 > 83 (J,Jy) peaks at around 
z = 0.6 — 1.0, and decays monotonically from z ~ 1 to 
z ~ 3. The shape of this decay is not reproduced by 
existing models of galaxy evolution. 

We have also obtained mid-infrared monochromatic 
and total infrared luminosities for all the sources, and 
have built luminosity functions at 12 fim. According to 
our results, the local luminosity function is relatively flat 
at faint luminosities (a ~ —1.2). Given the limitations in 
our data (in detection limit and areal coverage) , we esti- 
mated the luminosity functions in a number of ways that 
allowed us to understand the systematic errors. By fit- 
ting Schcchter (1976) functions and forms of the local lu- 
minosity function to the luminosity function data points 
in different redshift bins, we find: 1) the normalization of 
the luminosity function, <fr* , could be flat or increase at 
a maximum rate of (1 + z ) 2 - 1±0 - 6 up to z ~ 0.8; at higher 
redshifts, it seems to stay constant or decrease with red- 
shift; 2) the typical luminosity, L*, increases as at least 
(1 + Z y- &±1A to z ~ 0.8, and continues to increase at a 
roughly similar rate to higher z; 3) the best fits to our 
data predict an evolution where L* oc (1 + z) 3 0±0 - 3 and 
</>* oc (1 + z) 1 - 0±0 - 3 ; and 4) the low luminosity slope, a, 
is not well constrained between the values of a ~ —1.0 
to a ~ —1-7, but the best fits to our data indicate an 
almost flat value a < —1.3. 

We reproduce the previously seen rapid increase of the 
total infrared luminosity density of the Universe (and the 
cosmic star formation rate density) up to z ~ 1.4. This 
increase follows a (1 + z ) 4 0±0 - 2 law up to z — 0.8, with 
a declining rate up to z <~ 1.4. The slope at z < 1 is 
lower than that observed by some UV surveys, possibly 
indicating an evolution in the extinction properties of 
galaxies. At z > 1.4, we find no evidence for a decrease 
in the SFR density, but find a flat distribution up to z ~ 
3. Uncertainties in the faint end slope of the luminosity 
functions could affect these results significantly. 

Assuming an almost flat slope at faint luminosities 
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Fig. 11. — Relationship between the specific SFR (SFR per stellar mass unit) and the total stellar mass for the galaxies in our 24 fim 
survey. The sample is divided in 9 different redshift bins, each plotted in one panel. The lines, from left to right, correspond to constant 
SFRs of 1, 10, 100, and 1000 Mq yr _1 . Typical errors for each axis are also shown. 



for the luminosity functions at z > 0, our results in- 
dicate that the SFR density is dominated at low redshift 
(z < 0.5) by galaxies which are not very luminous in 
the infrared (£tir < lO 11 -^©)- The contribution from 
luminous infrared galaxies (10 11 < Ltir < 10 12 L Q ) in- 
creases rapidly from z ~ 0.4, forming approximately half 
of the total amount of newly-born stars by z ~ 0.7, while 
the starburst population declines steadily. At z = 1, Ul- 
traluminous Infrared Galaxies (Xtir > 10 12 L@) start to 
play a role, probably dominating the cosmic SFR density 
at z > 2. If we consider steeper values of the slope at 
faint luminosities for the luminosity functions at z > 0, 
the contribution to the total SFR density of starbursts is 
larger at high z, the evolution of LIRGs (relative to star- 
bursts) is not as marked as in the case of a flat slope, and 
all three galaxy types (starbursts, LIRGs, and ULIRGs) 
form approximately the same amount of stars at z ~ 2.5. 
The rapid increase of L* with z and our division of the 
cosmic star formation rate density according to the lumi- 
nosities of the contributing galaxies both agree. The role 
of ULIRGs in the overall star formation increases rapidly 
for z > 1.3. 



Finally, the distribution of masses and specific SFRs 
(SFR per stellar mass unit) of the galaxies in our survey 
seems to support a "downsizing" galaxy formation sce- 
nario, where the most massive galaxies would form first 
(z > 2), and the less massive systems would be continu- 
ously forming down to lower redshifts. 
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APPENDIX 

THE PHOTOMETRIC REDSHIFT TECHNIQUE 

Section 3 presented the main characteristics of our photometric redshift technique. Here we will describe the method 
in detail and discuss the quality of the redshift estimations. 

Data compilation 

Our photometric redshifts benefited from the use of a vast amount of data covering the UV, optical, NIR, and 
MIR spectral ranges. The main characteristics of each dataset, including the wavelengths, limiting magnitudes, and 
references for each filter, are given in Tables A3 and A4. 

Models and fitting technique 

We built two sets of SEDs from the 24 /im selected galaxies with spectroscopic data in CDFS (redshifts from VVDS) 
and HDFN (redshifts from TKRS and Cowie et al. 2004). In both cases, we only used the galaxies that were flagged 
as having accurate redshifts (confidence higher than 80%). We also selected only those galaxies with more than ten 
data points in their SEDs, assuring that the UV/optical, NIR, and MIR spectral ranges were covered. 

The observed SEDs were deredshifted (taken to z = 0) by transforming the effective wavelengths of the filters 
(calculated with the filter response and the detector quantum efficiency curves for each observation dataset 18 ) to 
the rest-frame and applying a (1 + z) factor to the flux density. The z = templates may suffer from inadequate 
K-correction. Indeed, when we convolved the filter curves with the z — templates and redshifted the results (to 
the original redshift of the galaxy), we did not recover the original observed SEDs. We used a minimization numeric 
method to obtain the best z = template compatible with the observed SED (i.e., we applied the K-correction using 
this minimization algorithm). 

We performed a visual inspection of all the templates, rejecting those with unusual shapes due to deviant photometry 
points. The final training sets were formed by 317 and 542 sources for the VVDS and HDFN datasets, respectively. To 
these, we also added the Devriendt et al. (1999) 17 empirically-calibrated templates. Some representative examples of 
our templates are shown in Figure A12. These examples show that the most prominent feature of the great majority 
of the SEDs is the 1.6 ^m stellar bump (marked with a dashed line). Nonetheless, for some galaxies this feature is 
hardly visible (see the bottom-right template in Figure A12). This absence of the stellar bump is observed in some of 

18 From this point, every time we refer to the filter response, we mean the filter transmission convolved with the detector response. 
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TABLE A3 

Characteristics of the data compiled for the CDFS. 
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NOTE. — (1) Name of the observing band. (2) Effective wavelength (in fim.) of the filter calculated by convolving the Vega spectrum (Colina & 
Bohlin 1994) with the transmission curve of the filtcr+dctcctor. (3) Limiting AB magnitudes defined as the third quartile of the magnitude distribution 
of our sample. (4) Source from where the data were obtained: a ESO Imaging Survey (EIS, Arnouts et al. 2002); b Las Campanas Infrared Survey 
(LCIS, Marzkc ct al. 1999); c The Great Observatories Origins Deep Survey (GOODS, Giavalisco ct al. 2004b); d EIS Deep Public Survey (EIS-DPS, 
Vandamc et al. 2001); e Classifying Objects by Medium-Band Observations -a spectrophotometric 17-filter survey- (COMB017, Wolf et al. 2004); f 
VIRMOS-VLT Deep Survey (VVDS, Lc Fcvrc et al. 2004). 



TABLE A4 

Characteristics of the data compiled for the HDFN. 
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wavelength (in /im) of the filter calculated by convolving the 
Vega spectrum (Colina & Bohlin 1994) with the transmission 
curve of the filter+detector. (3) Limiting AB magnitudes de- 
fined as the third quartile of the magnitude distribution of our 
sample. (4) Source from where the data were obtained: ™ pub- 
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the most infrared-luminous (presumably AGN-dominated) galaxies in the local Universe (see, e.g., Sanders et al. 1988; 
Devriendt et al. 1999). Despite the low resolution of the SEDs, some other spectral features are also visible, such as 
emission- lines in the UV/optical (marked with dotted lines), the 4000 A break (dashed-dotted line), or emission from 
PAHs in the MIR. 

After building the set of empirical templates, we redshifted them to values in the < z < 3 range using a step of 
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Fig. A12. — Example templates obtained from the VVDS dataset (data normalized to the 1 /im flux density). The original redshift of 
the source is given for each galaxy. Red stars show the original SED of the galaxy. Black open stars show the z = derived template 
without a K-correction calculation. Green stars show the recovered SED built by convolution of the z = non-K-corrected template with 
the filter response curves. Black filled stars (joined by a black line) show the final z = template after applying the K-correction. Blue 
stars show the final recovered redshifted template (which must coincide with the template formed by red stars). In each plot, the 1.6 /ivn 
bump is marked with a dashed line (and a magenta filled star), the Lyman and 4000 A breaks are marked with dashed-dotted lines, and 
the positions of the Lya, [Oil I] AA4959, 5007 and He? + [Nil] AA6548, 6584 emission-lines are marked with dotted lines. 



Az = 0.005. The redshift range was chosen based on expectations prior to launch (Dole et al. 2003). These redshifted 
templates were then convolved with the filter responses. 

The redshift estimation proceeded as follows. First, the entire observed SED was fitted with a Chebysev polynomial. 
When possible (i.e., when there were NIR and IRAC data and the stellar bump was present), the derivative of this 
polynomial was used to estimate the position of the 1.6 yum bump and its uncertainty. This position was used as a first 
guess and to constrain the final solution (in the range formed by the bump value and its error). We found this step 
to be a good procedure to get rid of outliers. Second, the observed and template fluxes were normalized to one of the 
bands (the reference band mentioned in Section 2.4) to account for the effects of different luminosities. Data points 
with dubious photometric calibration or repeated observations (two observations for the same filter) were removed 
before the fitting. In addition, the data points at the edges of the SEDs (the bluest and reddest filters, the latter 
always being 24 /im), were given smaller weights in the fitting, or removed when the templates had no data at those 
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Fig. A13. — Comparison between the photometric redshifts obtained in this work and the spectroscopic redshifts measured in CDFS (left 
panel) and in HDFN (right panel). Both comparisons refer to the photometric rcdshift results obtained with the complete template set 
(built with 317 sources from VVDS, 542 from TKRS and Cowie et al. 2004, and the 17 templates in Devriendt et al. 1999). Gray symbols 
are sources with unreliable spectroscopic redshifts. Open stars are sources detected in less than five bands. The dashed line shows the 
equality line, and the dash-dotted ones show the cr z /(l + z) < 0.1 area. 



wavelengths. Third, the templates and observed values were compared and a most probable redshift was calculated 
by minimizing a reduced % 2 estimator of the form: 

1 JVfiit (pi _ pi \2 

2 V ^ V template observed/ / . -t \ 

X ~ (A^ bscrV ed) 2 ( } 

where N^\ t is the number of filters considered, F t0 mpiate is the flux calculated for each redshifted template in the ith 
filter, and ^observed and Ai^ bserved are the measured fluxes and uncertainties in each filter. Errors in the redshift were 
calculated with a A\ 2 algorithm, and were quoted as the z-range for which the solution has a 68% probability of 
being correct. We only obtained redshifts for galaxies with more than four points in the SED (virtually all the sources 
mentioned in Section 2.4). 

Comparison with the spectroscopic sample 

The redshifts obtained with our photometric technique were compared with the spectroscopic values for all the 
galaxies in the sample. In Figure A13, we demonstrate that our procedure is able to recover the redshifts of the 
galaxies used as an input for the technique. 

The left panel in this Figure shows the comparison of photometric and spectroscopic redshifts for CDFS when using 
the complete template set (876 templates in total). The distribution of points is very symmetric around the equality 
line. The median value for the difference between the derived photometric redshift and the spectroscopic one (Sz 
hereafter) is Sz — 0.003, showing there is no systematic difference. Given that the wavelength scales in terms of (1 + z) 
when going to more distant sources, it seems more physically meaningful to discuss errors in terms of o~ z /(l + z) 
(Wolf et al. 2004), where u z is the absolute value of 8z. For the VVDS comparison, 90% of the objects have values of 
0-2/(1 + z) < 0.2, 82% of the objects have values of a z /(l + z) < 0.1, and 75% have cr z /(l + z) < 0.05. The average 
(median) a z j (1 + z) is 0.050 (0.012). These statistics improve if we only take into account sources with highly reliable 
redshifts (as stated by VVDS): 98% of the objects have values of a z j (1 + z) < 0.2, 92% of the objects have values of 
a z /(l + z) < 0.1, and 85% have a z /(l + z) < 0.05. 

The right panel of Figure A13 shows the comparison of our photometric redshifts with spectroscopic redshifts in the 
HDFN when using the complete template set. The statistics for this comparison (for the highly-reliable spectroscopic 
sample) are very similar to the CDFS case: < Sz > = —0.003, < a z /(l + z) >= 0.016, 98% of the objects have values of 
cr z /(l + z) < 0.2, 97% of the objects have values of o z j{\ + z) < 0.1, and 95% have a z j (1 + z) < 0.05. It is interesting 
to notice the smaller scatter of the points around the equality line in comparison with what we found in CDFS. This 
difference is due to the larger number of sources having NIR data in the HDFN, which makes the photometric rcdshift 
estimation more reliable because the 1.6 /xm bump position can be better constrained. 




photo photo 

Fig. A14. — Left: Comparison between the photometric rcdshifts obtained in this work and the spectroscopic redshifts for HDFN 
sources. This comparison refers to the photometric redshift results obtained by using the templates built in CDFS (317 templates). Right: 
Comparison between the photometric redshifts obtained in this work and the spectroscopic rcdshifts for CDFS sources. This comparison 
refers to the photometric redshift results obtained by using the templates built in HDFN (542 templates). For both panels, gray symbols 
arc sources with unreliable spectroscopic redshifts. Open stars are sources detected in less than five bands. The dashed line shows the 
equality line, and the dash-dotted ones show the u z /(l + z) < 0.1 area. 



The two panels in Figure A13 show some outliers, a total of 75 objects with a z j (1 + z) > 0.1 for the left panel (18% 
of the total sample of 425 sources in CDFS), and 22 objects with a z / (1 + z) > 0.1 for the right panel (4% of the total 
sample of 601 sources in HDFN). We will discuss these photometric redshift failures next, in order to characterize why 
the photometric redshift technique fails. None of these outliers were included within the template set, either because 
of a unreliable spectroscopic redshift, or because they did not have more than ten points in their SEDs. 

In CDFS, 60 out of the 75 outliers (i.e., 80% of the outliers) do not have highly-reliable spectroscopic redshifts (gray 
stars in Figure A13). Out of these, 26 present slight differences with the spectroscopic redshifts of <r z /(l + z) < 0.15 
(normally the redshift is overestimated). Out of the 34 remaining sources with unreliable spectroscopy, 18 present 
a clearly wrong photometric redshift. Half of these 18 sources present power-law like SEDs, and our method fails 
to obtain a good redshift. The other half do not present any particular problem, but the template selected by our 
technique gives a wrong redshift. The 16 remaining sources with non-reliable spectroscopy and cr z /(l + z) > 0.15 show 
SEDs clearly incompatible with the quoted spectroscopic redshifts. Out of the total 75 outliers in the left panel of 
Figure A13, 15 are secured spectroscopic identifications, with 10 of them presenting <r z /(l + z) < 0.15. The other 5 
do not have IRAC photometry, and our technique gives a very large redshift for them. 

In HDFN, 4 out of the 22 outliers have less than 5 points in their SEDs. This kind of object was removed from the 
photometric redshift sample in this paper due to the high uncertainties related to the small number of data points. In 
fact, there are a total of 7 galaxies within this group with reliable spectroscopic redshift in the entire HDFN sample, 
which gives more than 50% (4 out of 7) having a wrong photometric redshift. There are 8 more outliers with unreliable 
spectroscopic redshifts. Out of these, 3 present slightly overestimated photometric redshifts, with differences with the 
spectroscopic redshifts of <J z /(l + z) < 0.15. Another 2 sources (out of the 8 with unreliable spectroscopy) present a 
clearly wrong photometric redshift (one of the gray stars at z ~ 1.7, and the gray star at z p hoto — 0) due to highly 
deviant points in the SEDs (e.g., the F-band flux is 10 times larger than what would be expected based on other 
adjacent data points, the BvR bands, in the Zphoto ~ 1-7 case). This deviant point could be linked to source variability 
or deblending problems. The 3 remaining sources with unreliable spectroscopy show SEDs clearly incompatible with 
the quoted spectroscopic redshifts. Out of the other 10 outliers found in the right panel of Figure A13 (and having 
reliable spectroscopy), 4 of them present cr z /(l + z) < 0.15 (with a slightly overestimated photometric redshift). 
Another 6 are clearly wrong photometric rcdshifts, 3 of them lacking for NIR and/or IRAC data, and the other 3 
objects presenting very disturbed SEDs. 

In summary, the outliers in the two panels of Figure A13 are sources with wrong photometric redshifts due to 
photometry problems (about one third of them), objects with slightly overestimated redshifts (^50% of them), or 
probably wrong spectroscopic rcdshifts (10%-20%). This demonstrates that the fitting method does not introduce 
significant errors in the redshift determination. 

Most of the galaxies plotted in the left panel (more precisely, 60%) and right panel (90%) of Figure A13 were included 
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Fig. A15. — Comparison between the photometric redshifts obtained in this work (using the complete template set) and the ones obtained 
by COMB017 (Wolf ct al. 2004) in CDFS. Gray symbols are sources with Az > 0.05 as given by COMB017. Open stars are sources 
detected in less than five bands. The dashed line shows the equality line, and the dash-dotted ones show the cr z /(l + z) < 0.1 area. 



in the template set used to derive photometric rcdshift. Therefore, the numbers given are not entirely representative 
of the goodness of the method presented in this paper. A real test of the method is presented in Figure A14. Here we 
divided the sample in roughly two parts. One half of the sample was formed by the sources in CDFS, and the other 
half by sources in HDFN. We then obtained photometric redshifts for all the sources in one half of the sample by using 
the templates built from sources with a spectroscopic redshift within the other half of the sample. The left panel of 
Figure A14 shows the comparison of photometric and spectroscopic redshifts for the HDFN sources when using the 
templates built from CDFS sources. The right panel in Figure A14 shows the comparison for CDFS sources when 
only using templates from HDFN. The Figure demonstrates that we are able to obtain photometric redshifts (for the 
sources with secure spectroscopic redshift) with a z j (1 + z) < 0.1 for at least 80% of the sample (80% for the left panel 
and 84% for the right one), and redshifts with <r z /(l + z) < 0.2 for more than 90% of the galaxies (91% for the left 
panel and 94% for the right one). Other statistics are: < Sz > = —0.001, < cr z /(l + z) >= 0.078 for the left panel and 
< Sz >= -0.002, < a z /(l + z) >= 0.077 for the right panel. 

We also analyzed all the outliers on Figure A14 to characterize why the photometric redshift technique fails. In the 
left panel, there are 601 galaxies with spectroscopic redshifts, 125 (21% of the total spectroscopic sample in HDFN) 
of which present <r z /(l + z) > 0.1. Out of this number, 67 sources (54% of the outliers) have <J Z /(1 + z) < 0.2, half 
of them with overestimated photometric redshifts, half with underestimated values. Within the 58 sources remaining 
(from the 125 outliers), 7 galaxies present very disturbed SEDs (6% of the outliers). Another 36 objects are clear 
photometric rcdshift errors (29% of the outliers), probably because of the lack of NIR and/or IRAC data (13 objects 
out of the 36), or because they present power-law SEDs (10 objects). The rest, 15 sources (12% of the outliers), are 
probable spectroscopic redshift errors (virtually all of them flagged as unreliable spectroscopic estimations). In the 
right panel of Figure A14, there are 425 galaxies, 81 (19% of the total spectroscopic sample in CDFS) which present 
er z /(l + z) > 0.1. Out of this number, 25 (31% of the outliers) have cr z /(l + z) > 0.2, and only 10 of those are 
labeled as secure spectroscopic redshifts, all with very disturbed SEDs. In summary, less than 20% of the galaxies 
present <r z /(l + z) > 0.1, at least one third of those present (T z /(1 + z) > 0.2, and from these 10%, at least one 
fourth are probable spectroscopic redshift errors. Most of the sources with a z /(l + z) > 0.1 either lack for NIR/IRAC 
photometry, or present disturbed SEDs probably linked to source variability (related to AGN activity) or dcblending 
problems. In fact, ~6% of the sources with a z /(l + z) < 0.1 present multiple identifications (within the search radius), 
in comparison with a ~10% for galaxies with c z /(l + z) > 0.1. 

Finally, we also compared our results with other photometric redshift surveys such as COMB017 (Wolf et al. 
2004). COMBO 17 and our redshifts agree very well up to z ~ 1, with more than 90% of the sources being within 
the <7 Z /(1 + z) < 0.1 area (just for a highly-reliable sample, i.e., sources with Az < 0.05 as given by COMB017). 
However, there are two issues in the comparison, that we can clearly see in Figure A15. First, the distribution of 
points is not symmetric at z > 1. Many sources lying near the COMB017 redshift limit (z ~ 1.4) are placed at higher 
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rcdshifts by our photometric redshift method, most of them presenting high COMBO 17 redshift uncertainties. For 
< zcoMBon < 1-4, we find: < Sz >= 0.052, < a z /(l + z) >= 0.114, and 75% of the sources have a z /(l + z) < 0.1. 
COMB017 is known to have a deficiency when estimating redshifts at z > 1, because the useful spectral features go 
out of their optical filter set. These are also the faintest sources in the COMB017 sample, and the photometric redshift 
method gets more uncertain. In our case, NIR data are available for 50% of the sources, and IRAC photometry for 
virtually all of them. These bands allow us to trace the spectral features used by COMB017 to high redshift (z ~ 3), 
and to break redshift degeneracies coming from the misidentification of the Balmer break with the Lyman break. The 
second issue seen in Figure A15 is that some galaxies (approximately 5% of the common sources between our survey 
and COMB017) arc placed at z < 0.2 by COMB017 and at higher redshifts by our work. Le Floc'h et al. (2005) 
also noticed this effect in studying the luminosities of MIPS sources in CDFS, concluding that they must be redshift 
outliers in the COMB017 survey. 

The statistics and the fraction of outliers of our photometric redshift technique are comparable to most other 
photometric redshift works in the literature (e.g., Connolly et al. 1995; Chen et al. 2003; Firth et al. 2003; Rowan- 
Robinson 2003; Babbedge et al. 2004; Zheng et al. 2004), and only slightly worse than some surveys aimed at obtaining 
high-quality photometric redshift by using adequate sets of narrow-band filters, although the new Spitzer data helps 
to obtain more reliable and accurate redshifts for z > 1 sources (e.g., Wolf et al. 2004). 

Although Le Floc'h et al. (2004) demonstrated the ability of Spitzer to identify sources in the 'redshift desert' 
(1.5 < z < 3.0, Steidel et al. 2004), our photo-z technique might have deficiencies for this redshift range. The 
spectroscopic surveys are highly biased towards z < 1 galaxies, given their technical limitations for very faint objects. 
Indeed, our template set does not have many galaxies lying at high redshifts: 80% of the templates correspond to 
z < 1 sources, 16% to 1.0 < z < 1.4 galaxies, and only 4% above z = 1.4. Our technique assumes that the most distant 
galaxies in our sample are very similar (in their SED properties) to the closer ones, i.e., we can find templates fitting 
the z > 1.4 galaxies among the z < 1.4 sample. Figure 5 shows that this extrapolation is not unreasonable, given 
that the LIRG population (probably showing similar properties at all redshifts) dominates our sample (in number of 
sources) from z ~ 0.8 up to z ~ 2.0. Only above z ~ 2.0, the ULIRGs are statistically relevant in number, and we do 
not have many templates with known spectroscopic redshifts. 

Having this possible problem in mind, we performed a visual inspection of the fits for the galaxies identified as being 
at z > 1.4 by our photo-z technique. This test is represented in Figure A16, which plots 10 sources selected randomly 
from the entire sample in the redshift range < z < 1.4 (left column), and 10 more at 1.4 < z < 3.0 (right column). 
This figure gives an overview of the reliability of our photo-z technique. Most of the galaxies in our sample show a clear 
1.6 /im stellar bump, a feature that allows a reliable photometric redshift determination. Indeed, all z < 1.4 galaxies in 
Figure A16 show a negative average slope in the IRAC bands, and a positive slope in the optical/NIR bands, pointing 
to a clear 1.6 /<m bump. The determination of the exact position of this bump is, however, importantly affected by the 
availability of NIR data (and consequently, the uncertainty in the photometric redshift). For the sources at z > 1.4, at 
least 7 galaxies present a change in slope (from positive to negative) inside the spectral range covered by IRAC, which 
indicates the presence of the 1.6 ^m bump (for a 1.3 < z < 4 galaxy). The visual inspection of the randomly selected 
sources in Figure A16 revealed that the assigned redshift was dubious for 20%-25% of the sources, all of them lacking 
a marked 1.6 /im stellar bump (probably linked to the presence of an AGN), and preferentially lying at z > 1.4. This 
percentage is very similar to the total reliability of our technique (<~ 80%), previously discussed in this Section, which 
suggests that the procedure remains applicable to the highest redshift range. A better coverage of the 'redshift desert' 
for LIRGs and ULIRGs by spectroscopic surveys would be desirable to reduce the uncertainty of the results achieved 
at z > 1.4. 

SYSTEMATIC UNCERTAINTIES IN THE LUMINOSITY FUNCTION ESTIMATION LINKED TO 

PHOTOMETRIC REDSHIFT ERRORS 

In this Appendix, we will investigate the effect of the photometric redshift errors discussed in Appendix A (i.e., 
scatter and fraction of outliers) in the estimation of the luminosity functions. 

The photometric redshift errors propagate in the calculation of luminosities. Given that the estimation of the 
luminosity functions involves binning and averaging of luminosities, the uncertainties linked to redshift errors could 
in principle diminish, and they should most probably be random. However, there could also exist some systematic 
uncertainties, specially in the extremes of the luminosity function (faint and bright ends), where the number of detected 
galaxies is small, and the redshift outliers are preferentially found. 

To understand the uncertainties in the luminosity function linked to photometric redshift errors, we performed a 
Monte Carlo simulation similar to the one presented in Chen et al. (2003). This simulation consisted in the genera- 
tion of an artificial catalog of galaxies following an assumed input Schechter luminosity function and presenting the 
characteristic limiting fluxes of our survey. Each galaxy was assigned a random redshift between z = and z = 3. 
This redshift was perturbed by an amount linked to the probability distribution built from the scatter of points in the 
comparison between spectroscopic and photometric redshifts presented in Figure A14. This probability distribution 
of photometric redshift uncertainties accounts for the scatter of points around the equality line in Figure A14, and for 
the outliers (i.e., the most deviant points in that scattering plot, caused by photometry problems or any other issue). 
Finally, we estimated the luminosity function for three representative redshift intervals: 0.0 < z < 0.2, 0.8 < z < 1.0, 
and 1.8 < z < 2.2. We will concentrate our discussion in the lowest redshift interval, 0.0 < z < 0.2, where we can 
explore the largest luminosity range. Note that this interval only probes the photometric redshift errors toward higher 
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Fig. A16. — Some randomly-selected examples of the fits obtained in the photometric redshift estimation. On the left, we show ten 
sources selected randomly in the redshift range < z < 1.4, and on the right ten more randomly selected sources at 1.4 < z < 3.0 are 
given. Photometry points are plotted with stars, and the best template fit to the data is shown with a solid line. 



rcdshifts. Given that we obtained very similar results for the other intervals, we concluded that the errors toward 
lower redshifts do not affect the results from our simulation significantly. 

The results from the Monte Carlo simulation are plotted in Figure B17. The input luminosity function (continuous 
line) is not well recovered with the standard SWML method (open stars fitted by the dotted line). As we previously 
suggested, the photometric redshift errors affect the faint and bright ends of the luminosity function, resulting on an 
overestimation of the galaxy density in these ranges, which turns into an overestimation of both a (by ~20%) and 
L* (by ^0.3 dex). This result is consistent with that found by Chen et al. (2003). Note, however, that although the 
luminosity function parameters are not properly recovered due to photometric redshift errors, the integrated luminosity 
density does not change significantly (less than 0.1 dex, below the estimated uncertainty 0.25 dex). 
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Fig. B17. — Results from the Monte Carlo simulation performed to investigate the systematic uncertainties in the 12 fim luminosity 
function estimation introduced by the photometric redshift errors. The solid line shows the input Schechter luminosity function, whose a 
and L* values are indicated by the triangle in the inset. The open stars and dotted line (also in the inset, with the contour delimiting 
the 99.99% probability area) show the recovered luminosity function for the artificial catalog. These open stars have been offset to higher 
luminosities by a constant small amount for clarity. The filled stars and dash-dotted line (also in the inset) represent the final luminosity 
function obtained with the modified SWML method (using a Monte Carlo iterative method). This technique accounts for the typical 
luminosity function errors (linked to the number of galaxies in each luminosity bin), as well as for redshift uncertainties. 



To cope with the previously described systematic errors, we modified the SWML luminosity function method to 
include the photometric redshift and luminosity uncertainties. We used a Monte Carlo approach where each redshift 
was treated as statistical variable with an average and an uncertainty. The average was the photometric redshift 
given by our technique and the uncertainty was derived from the probability distribution of redshift errors extracted 
from Figure A14. The redshifts for the whole catalog of galaxies in our survey were randomly perturbed according to 
this distribution, and the 12 fim luminosities and their uncertainties were recalculated. The luminosity function was 
estimated using the SWML method for the new catalog. We iterated this process 1000 times, and the final luminosity 
function points (given in Figures 6 and 7), the fitting parameters, and the associated errors for all these quantities 
were derived from the distribution of solutions. 

The results from the modified SWML method for the simulation described previously are shown in Figure B17 with 
filled stars and a dash-dotted line (also in the inset). We are able to recover the input luminosity function parameters 
with high accuracy, reducing the systematic errors to less than 2% for a, and less than 0.04 dex for both L* and <f)* . In 
general, the modified SWML method obtained luminosity functions with smaller faint-end slopes and lower L* values 
than the traditional SWML technique for all the redshift ranges. 



